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Abstract 

The cortex is composed of large-scale cell assemblies sharing the same individual properties 
and receiving the same input, in charge of certain functions, and subject to noise. Such 
assemblies are characterized by specific space locations and space-dependent delayed inter- 
actions. The mean-field equations for such systems were rigorously derived in a recent paper 
for general models, under mild assumptions on the network, using probabilistic methods. 
We summarize and investigate general implications of this result. We then address the dy- 
namics of these stochastic neural field equations in the case of firing-rate neurons. This is 
a unique case where the very complex stochastic mean-field equations exactly reduce to a 
set of delayed differential or integro- differential equations on the two first moments of the 
solutions, this reduction being possible due to the Gaussian nature of the solutions. The 
obtained equations differ from more customary approaches in that it incorporates intrinsic 
noise levels nonlinearly and make explicit the interaction between the mean activity and 
its correlations. We analyze the dynamics of these equations, with a particular focus on 
the influence of noise levels on shaping the collective response of neural assemblies and 
brain states. Cascades of Hopf bifurcations are observed as a function of noise amplitude, 
for noise levels small enough, and delays, in a finite-population system. The presence of 
spatially homogeneous solutions in law is discussed in different non-delayed neural fields 
and an instability, as noise amplitude is varied, of the homogeneous state, is found. In 
these regimes, very complex irregular and structured spatio-temporal patterns of activity 
are exhibited including in particular wave or bump splitting. 

Keywords: mean-field equations, neural fields, bifurcations, integro-differential equations, 
noise, infinite-dimensional dynamical systems, propagation of chaos, delayed stochastic 
integro-differential equations. 



Contents 

1 Introduction 2 

Email address: jonathan.touboul(9college-de-f rance.fr (Jonathan Touboul) 

Preprint submitted to Elsevier November 1, 2011 



2 Mathematical results on mean-fields models for neural fields 4 

2.1 Finite number of populations 4 

2.2 Infinite number of populations 6 

3 Mean- field equations for firing-rate models 12 

3.1 Finite-population Wilson and Cowan Equations 13 

3.2 Spatially extended neural fields 16 

4 Analysis of the solutions 21 

4.1 Finite-populations networks 21 

4.2 Neural Fields Dynamics 25 

4.2.1 A single layer neural field: Turing-Hopf instabilities and noise 26 

4.2.2 Dynamic Turing patterns in a two layers neural field with periodic 
boundary conditions 30 

5 Discussion 40 

Appendix A Existence and Uniqueness of solutions of the synchronized 

mean-field solution 43 

Appendix B Spatio-temporal patterns for one-dimensional neural fields 

with reflective or zero boundary conditions 45 

Appendix B.l Reflective boundary conditions 45 

Appendix B.2 Zero boundary conditions 45 



1. Introduction 

The activity of the brain is often characterized by large-scale macroscopic states result- 
ing of the structured interaction of a very large number of neurons. This interaction yield 
meaningful signals accessible from non-invasive imaging techniques (EEG/MEG) and often 
used for diagnosis. Finer analysis of the brain's constitution identifies anatomical struc- 
tures, such as the cortical columns, composed of the order of few thousands to one hundred 
thousand neurons belonging to a few different species, in charge of specific functions, shar- 
ing the same input and strongly interconnected, communicating through the emission of 
action potentials after delays depending on the anatomical or functional distance between 
themselves. Two interesting well-understood cortical areas are the primary visual cortex 
of certain mammals and the rat's barrel cortex. In the primary visual cortex VI, neurons 
can be divided into orientation preference columns responding to specific orientations in 
visual stimuli, forming specific patchy connections [1, 2]. Similarly, the rat's barrel cortex 
presents a clear columnar organization, and neurons responding to the sensory information 
of a particular whisker anatomically gather into barrels [3, 4]. More generally, it is now 
understood that brain's activity results from the interactions of different cells, among which 
neurons, manifesting highly complex behaviors often characterized by the intense presence 
of noise. The communication between two such columns is characterized by a delay due to 
the transport of information through axons at a finite speed and to the typical time the 
synaptic transmission takes. These delays have a clear role in shaping the neuronal activity, 
as established by different authors (see e.g. [5, 6, 7, 8]). Several relevant brain states rely 
on the coordinated behaviors of large neural assemblies, and recently raised the interest of 
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physiologists and computational neuroscientists, among which we shall cite the rapid com- 
plex answers to specific stimuli [9], decorrelated activity [10, 11], large scale oscillations [12], 
synchronization [13, 14], and spatio-temporal pattern formation [15, 16, 17]. 

Most computational models of neural area have relied on continuum limits ever since 
the seminal work of Wilson, Cowan and Amari [18, 19, 20, 21], in which the activity is 
represented through a macroscopic variable, the population-averaged firing rate, that is 
generally assumed to be deterministic. This approach successfully accounted for several 
phenomena, for instance for the problem of spatio-temporal pattern formation in spatially 
extended models (see e.g. [16, 22, 23, 24]). In this approach, the effect of noise is neglected 
in large populations. Increasingly many models now consider that the different intrinsic 
or extrinsic noise sources are a meaningful part of the neuronal signal rather than a pure 
disturbing effect, that conveys important information [25]. For instance, sparsely connected 
neural networks were studied, in wihch the sparse connectivity assumption allows studying 
regimes where the activity is uncorrelated [26, 27, 28]. Again, the emergent global activity 
of the population in the limit of an infinite number of neurons is deterministic, and evolves 
according to a mean- field firing rate equation. The stochastic nature of the firing was also 
developed through such methods as the population density [29] , or from a Markovian model 
governing the firing dynamics of the neurons in the network, where the transition probability 
satisfies master equation. This modeling recently gathered the interest of the community (see 
e.g. [30, 31, 32, 33, 34, 35]). Most of these approaches are proved correct in some parameter 
regions using statistical physics tools such as path integrals and Van-Kampen expansions, 
and involve a moment expansion and truncation. 

We are interested here to start include noise at the level of each individual neuron. From 
the mathematical viewpoint, each neuron is a diffusion process, and different neurons belong 
to particular populations in which they receive a common input and are interconnected in 
a similar fashion. These equations motivated by biological analysis, are studied in a prob- 
abilistic setting in a recent paper [36], taking into account spatial information and delayed 
interconnections, using probabilistic methods from the propagation of chaos domain (see 
e.g. [37]). The results of this analysis are summarized in section 2. Under mild assumptions 
on the neuronal dynamics, the regularity of the interactions and the number of neurons per 
population, this analysis provides a mean-field equation corresponding to the behavior of a 
particular neuron in the network. These equations, in the original general framework, apply 
to such neurons as the Hodgkin-Huxley and Fitzhugh-Nagumo neuron models, which are 
relevant representations of individual neuronal activity. We analyze here the implications 
of this result in that general setting in the discussion section. We also discuss an extended 
notion of synchronization, the synchronization in law^ and provide general conditions for 
this property to occur in spatially extended mean-field neural fields systems. The mean- 
field equations obtained in [36] are implicit equations on probability distributions, and as 
such extremely hard to tame in general. Simulations of such equatiosn are relatively hard 
to perform, as illustrated in the case where no delay and no spatial extension is taken into 
account [38], and this difficulty makes it hard to identify qualitative effects of noise levels in 
such systems. In order to start uncovering the structure of the solutions of these equations 
and possible effects of noise levels each individual neuron is submitted to, we consider in 
section 3 a simplified neuron model based on the description of firing-rates. In this par- 
ticular case, we demonstrate that the solution of the complicated mean- field equations are 
Gaussian, similarly to what was done in a way simpler model in [39]. Characterizing the 
mean and covariance functions are hence sufficient to describe the dynamics of the solu- 
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tions, and these are shown to exactly reduce to a set of deterministic delayed differential 
or integro-differential equations. The derivation of these equations and the bifurcations of 
these equations are studied, in the finite-population case (section 3.1) and an intricate bi- 
furcation structure with the appearance of a cascade of Hopf bifurcations as a function of 
noise levels and delays, that tend to accumulate as delays are increased in a restricted region 
of noise levels. This bifurcation structure results in the presence of an increasing number 
of unstable cycles that provoke extremely complex and irregular transient behaviors. In 
the continuous-populations case (section 3.2), we show noise levels very strongly shape the 
form of the spatio-temporal solutions of the system. Spatially homogeneous stationary or 
periodic solutions are found as a function of the noise levels. An additional state is exhibited 
for small values of the noise levels, in which spatially inhomogeneous, non-stationary spatio- 
temporal patterns with a complicated shape arise. These qualitative results are reproduced 
for different unidimensional neural fields with different topologies and connectivities, and 
either satisfying or not the synchronization in law condition, in appendix Appendix B. 
Implications of these results are eventually discussed in section 5 together with some open 
problems the present study motivates. 

2. Mathematical results on mean-fields models for neural fields 

In [36] a rigorous derivation of the mean-field equations is performed for a wide class of 
models encompassing classical neuron models such as the Hodgkin-Huxley or the Fitzhugh- 
Nagumo models, taking into account neural fields spatial structures and space-dependent de- 
lays. The general results obtained in that article are summarized here. In all the manuscript, 
we are working in a complete probability space (1], J^, P) satisfying the usual conditions, and 
endowed with a filtration (^t)^- The state of each neuron i in the network is described by 

a d-dimensional variable X* G £^ =^R^, typically corresponding to the membrane potential 
of the neuron and possibly additional variables such as those related to ionic concentrations 
and channels. We distinguish the cases where the number of populations is finite and the 
case where it is infinite corresponding to the spatially extended neural fields. 

2.1. Finite number of populations 

We start by summarizing the results of [36] in the case where the network is composed 
of a finite number of populations a G {1, • ' ' ^P}- Let p denote the population function that 
associates to a neuron i the population a it belongs to: p{i) = a, and Na{N) the number 
of neuron in population a for a network size N (we will generally drop the dependence of 
Na in N since no confusion is possible), and assume that this size tends to infinity in the 
limit N ^ oo. The state of a neuron i in population a, X* is characterized by its intrinsic 
dynamics governed by a drift function Ga : x E E (including the intrinsic dynamics 
and the deterministic inputs) and a diffusion matrix - x E R^. 



where (ly/)^^^ constitutes a sequence of m x d-dimensional adapted Brownian motions. 
Neuron j of population 7 produces an infinitesimal current received at time t by the neuron 
i modeled as: 



dXi = Ga{t, Xi) dt + gc,{t, Xi) dW^ 



t 5 
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The terms rja-f and /la-f denote signed finite measures modeling the delays occurring in the 
information transmission through the axons, and /^q,^ are smooth functions of ExE E 
and (5j'^)iGN,7=i, -- ,p are d x d independent adapted Brownian motions. In this setting, 
the network equations read: 



p N 



dXl^"" = G^{t,Xl^'')dt^Y. W I ^^^i^t''^^ifs)dVa,{s)dt 

P . N 

+ g^{t,Xl^'^)-dWi^^— ^ / Pa.iXl^'^.Xl 
7=1 ^ .7 = 1,«(.7)=7'^"^ 



P , N 

'•^ ^/;^)dMa7(«)-rf5r- (1) 



i=i,p(i)=7 ' 

For the sake of simplicity, we will denote in the sequel the sum over neurons of population 7, 
XljLi p(j)='y7 simply by X^j^Zi- These equations are well-posed stochastic differential equa- 
tions on the infinite-dimensional space Cr = r, 0],£^) (see e.g. [40, 41]) of continuous 
functions of [— r, 0] with values in E. We will generally make the assumption that the net- 
work has chaotic initial condition, i.e. independent initial conditions, identically distributed 
for neurons belonging to the same population. In details, for {(Q{t),a = 1 • • • P) G C^- a 
stochastic process with independent components, a chaotic initial condition on the network 
consists in setting the initial condition of all neurons i of population a to an independent 
copy of Co • The following mild technical assumptions on the functions governing the net- 
work's dynamics: 

(HI). Ga and ga are uniformly locally Lipschitz-continuous with respect to their second 
variable 

(H2). and P^-f are L-Lipschitz-continuous, i.e. there exists a positive constant L such 
that for all (x, y) and {x' ^y') in E x E we have: 

|e«^(x, y) - e«^(x^ ^01 <L{\x-x'\^\y- y'\) 

for 6 G {6, (3}. 
(as). There exists a ^ > such that: 

max(|6,^(x,z)|M/3,^(x,z)p) < ^(1 + |xp) 

(114). The drift and diffusion functions satisfy the monotone growth condition: 

X^fa{t,x) + ^\ga{t,x)\^ <K{1 + 

Under these assumptions, the following theorem is proved in [36]: 

Theorem 1. Let us consider a network composed of N neurons belonging to P popula- 
tions. Assuming the regularity assumptions (H1)-(H4) are satisfied and that each neu- 
ron of the network chaotic initial conditions (with law (""'^ for neurons belonging to the 
same population a). Fix I G N* and {ii,-- - ,ii} I neurons of the network. The law of 
(X^^'^,--- r < t < T) solution of the network equations (1) converges towards 

jyipM (g) . . . (g) jyipM when N ^ 00 where is the law of the unique solution X of the 
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equation: 

dX^=Go.{t,X^)dt^Y. / ^AK^{X^X%s)]dria^{s)dt 

^go,(t,X^)dW^^Y. ^APo.p{X^Xis)Wo.p[s)dBr (2) 

with initial condition (Cf,Q^ = 1---P, t G [— t, 0]). In equation {2), the process (Y^) is a 
process independent of, and has the same law as and Ey denotes the expectation with 
respect to the process Y. In other terms, if ml denotes the probability distribution of X] , 
the mean-field equations can be equivalently written as: 

dX^=Go.{t,X^)dt^Y. f I K^{K.y)rnlUdy)dria^{s)dt 

-^g^{t,X^)dW^^Y. f I Pa^{X^.y)m]^,{dy)dfi^^{s)dBr- 

2.2. Infinite number of populations 

We now summarize results of [36] in the case of neural fields spanning over an infinite, 
continuous compact space F, modeling a spatial or functional extension of a cortical struc- 
ture. These structures are made of several sub-areas, the neural populations, composed of a 
large number of strongly interconnected neurons. The set T considered is a finite-dimensional 
compact set. Typically, it is considered, when populations are described by their location 
on the cortex, that F as a compact subset of for some q G N* (in Appendix B, we will 
consider in particular the case F = [0,1]). When considering that populations are defined 
by the neuron's function, the shape of F can take different forms depending on the geometry 
of the feature space. For instance in the case of the modeling of pinwheels primary visual 
area, neurons are indexed by their preferred orientation that can be represented in the torus 
F = (this is the case of section 4.2). 

In that case, the total number of populations is a diverging function P{N) of the total 
number of neurons and populations densely cover F in the limit N ^ oo with a certain 
density. The populations are characterized by their location (ri,--- ,rp(7v)) ^ F^*^^^ (see 
figure Fig. 1), and by the number of neurons in each population in a network of size TV, 
{Ni{N),...,Np(^N^{N)) (we hence have Yl^^^ ^-fi^) = Populations are randomly 
and independently drawn in a law A(-)/A(F) where A is a square integrable^, finite density 

over F (i.e. a positive function such that Jp A^(r) dr =^(A^)(F) < oo), accounting for possible 
inhomogeneities of the neural field. Note that the square integrability condition readily 
implies that Jp A(r) =^ A(F) < oo. This degree of freedom allowed for accounting for 
different topologies and geometries of the physical space or feature space (for instance, a 
recent model of texture neural field involves a feature space having an hyperbolic geometry. 



-"^Only integrability was assumed in [36], the square integrability will be needed for our developments on 
firing-rate models. 
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Figure 1: A typical architecture of neural field: cylinders represent neural populations as cortical columns 
spanning across the cortex. Neuron j in the green population at location rj G F communicates with neuron 
i in the red population at location r^, and the information sent is received with a delay r(ri,rj) depending 
on the spatial location of each neuron, and with a synaptic weight J(ri,rj) also depending on the neural 
populations of i and j. 



see e.g. [42, 43]). The expectation over the reahzations of the space locations {ra) is denoted 



As observed in [36], there is a critical competition between the number of populations 
and the total number of neurons. Indeed, in order for averaging effects to occur in the 
neural field, a very large number of neurons in each population is required, competing with 
the number of populations, which is required to fill the space P. A technical assumption 
made in [36] was the following: 



This assumption ensures heuristically most populations are made of a diverging number of 
neurons. This interpretation is relevant for our neuronal populations problem since, it is un- 
derstood that populations roughly contain the same number of neurons and that this number 
of neurons is orders of magnitude larger than the number of populations (see e.g. [44]). In 
order to comply with more standard approaches in neural field theory, the functions Ga{t, x) 
are replaced by G{ra,t,x) and similarly ga{t^x) by g{ra^t,x). The number of populations 
diverging, it is necessary in our description to scale the synaptic weight as the number of 
populations increases, in order to ensure the finiteness of the input to one single neuron^, and 
define the interaction functions between population a and population 7 at location r^, de- 
noted in the finite-population case ba,^{x,y) (resp. Pa,-f{x^y)) by p(V) ^(^q^ v) (I'^sp. 
p(Ar) ^(^Q;^ r^, X, ?/)). The delay measures r]{ra^r^^s) and jiira-ir^^ s) typically involve the 
distance Htc^ — r^||r and the transmission velocity of the information in the axons. These are 
again assumed to be finite signed measures, charging the interval [— r, 0], with total variation 
uniformly (in Tq, and r^) bounded by a quantity denoted n. 

We define a sequences of independent Brownian motions (W/) in A^(C([0, T]), i^^x^) and 
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(3) 



^This is due in particular to the limited resources of the neuronal environment. 
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a sequence of independent cylindrical Brownian motions Bl{r') in A1(C([0, T]), ]L^(r, R"^^^)) 
for i G N, i.e. centered Gaussian processes with covariance function E [Bl(r)Bi(r')\ =tAs 
if i = j and r = r' ^ and otherwise (see [40] for existence and properties of this process). 
The stochastic network equations in the read, for p{i) = a and a network of size N: 



^ ' 7=1 j=l T -^-T 

PiN) N, 

+ g{r^,t,Xl''')-dWi + -^ (i{ra,r^,Xl^'',Xlf,)dii{r^,r^,s)-dBi{r^). 

■^V^^J j.^l ^^7 J-T 

(4) 



They make the following suitable regularity assumptions on the functions governing the 
dynamics of the neurons: 

(GHl). G(r, •) and ^(r, •) are uniformly locally Lipschitz-continuous, 

(GH2). 6(r, rVr) and /3(r, rVr) are uniformly Lipschitz-continuous, i.e. there exists a pos- 
itive constant L such that for all (x^y) and {x'^y') in E x E and all (r, r') G we 
have: 

|6(r, X, y) - e(r, x^ y') \ < L {\x - x'\ ^ \y - y'\) 

for e e {6, 
(GH3). There exists a > such that: 

max(|6(r, r', x, z)f, |/?(r, r', x, z)f) < K{1 + 

(GH4). The drift and diffusion functions satisfy, uniformly in space (r) and time (t), the 
inequality: 

x^/(r,t,x) + i|^(r,t,x)p < (1 + 
Under these assumptions, it is proved in [36] the following: 

Theorem 2. Let us consider a network composed of N neurons belonging to P{N) popu- 
lations such that the assumptions {?>), (GH1)-(GH4) are satisfied and that each neuron of 
the network has independent initial conditions, that are identically distributed population per 
population (with law C^{r) for neurons belonging to the population located at r G T). Fix 
I e N* and {n, • • • ^ii} I neurons of the network. The law of {Xl^' , • ' " ? -^t'' t < T) 

solution of the network equations (4) converges towards mt{rp(^i^^) • • • (g) mt{rp(^i^^) when 
N oo, where mtir) is the law of Xtir), the unique solution of the mean-field equations: 

dXt{r) = G{r,t,Xt{r))dt+ [ [ lEz[b{r,r' ,Xt{r), Zt+s{r'))]d7]{r,r' , s) X{r')dr' dt 



+ g{r,t,Xt{r))-dWt{r)+ J E2[/^(r, r', Xt(r), Zt+,(r'))]dM'', r', s) • dBt(r, r')A(r') rfr' 

(5) 

with initial condition (C{^(r),r G r,t G [— t, 0])^ where the process {Zt{r)) is independent of 
and has the same law as Xt{r). 
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These results presented in the summary can appear very formal due to the generality of 
the models considered and to the mathematical approach developed of [36]. However, even at 
this level of generality, one can deduce several implications of these properties with regards 
to neuroscience applications (see section 5). It is important to realize at this point that 
even if the mean-field approach reduces an infinite system of interacting diffusion processes 
into a single well-posed equation, a major issue one faces is the concrete identification, 
characterization and simulation of the solution of the mean-field equations. In our spatial 
and delayed setting, this concern is even more true. We show in this paper that there exists a 
particular model of neuronal dynamics related to the so-called firing-rate neurons, for which 
mean-field equations are rigorously reducible to a set of deterministic equations that are 
amenable to analytic treatment (see sections 3 and 4). We emphasize the fact that firing- 
rate neurons, though highly popular and widely used in the computational neuroscience 
community, are not the most relevant from the biological viewpoint, and this restriction 
justifies the necessity of dealing with such complex models as those described in [36]. 

A phenomenon of great interest in neuroscience is the synchronization, understood as the 
fact that neurons manifest the same behavior after a transient phase, or the polychronization 
corresponding to the fact that the neural field form a few synchronized clusters (see e.g. [14]). 
Our stochastic setting suggests to extend this notion to the wider notion of synchronization 
in law, corresponding to the fact that the law of state of different neurons or different 
populations converge (as time increases) towards the same probability distribution. This 
notion of synchronization generalizes the notion of synchronized oscillations to solutions 
that are stationary in time, or that present specific time evolutions. In particular, spatially 
homogeneous in law solutions of the neural field equations will be said synchronized in law. 
It is clear from the results of theorem 2 that different neurons of the same population in the 
network have, in the mean-field limit, the same probability distribution provided that they 
had iid initial conditions, hence synchronize in law, and that different neurons belonging to a 
few different populations polychronize in law, forming clusters depending on the population 
the neurons belong to. A more complex question is the synchronization or polychronization 
in law between different populations. The following proposition provides sufficient conditions 
for solutions to be spatially homogeneous in law in the neural field setting. To this end, we 
extend the simple sufficient conditions of [45] in the stochastic mean-field setting. We have 
the following: 

Proposition 3. Assume that the distribution of the initial condition C?(r) is chaotic and 
does not depend on r eT, and that the functions G{r, t, x) and g{r, t, x) does not depend on 
r. Assume moreover that the law of the quantities: 



do not depend on r for any ip a measurable function (this ensures that the global input a 
neuron receives is independent of its spatial location). Then the solution of the mean-field 
equation (5) is spatially homogeneous (i.e. does not depend on r, or is fully synchronized). 
Moreover, if we fix ro G V, then for each r ^T, the law of Xt{r) solution of the mean-field 




(6) 
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equations is equal to the law of the solution of the equation: 

Xt{ro) = Xo{ro) ^ ds(^G{ro,s,Xs{ro)) ^l^z[B{ro,^ 

+ Ez[i/(ro,X,(ro),Z(.)(ro))]+ f dWs{r)g{ro^ Xs{ro)). (7) 

Jo 

which has a unique solution. 

Remark. Under the spatial homogeneity condition of proposition 3, two neurons belonging to 
two different populations r and r' are heuristically interchangeable in the limit N ^ oo, and hence 
the finite-population setting heuristically applies with only one population network, made of neurons 
with intrinsic dynamics function G{ro,s,x) and diffusion function g{ro, s, x) , deterministic (resp. 
stochastic) interaction functions i?(ro, x, (ro)) (resp.H(ro,x, Z(^.^{ro))). 

Proof. The existence and uniqueness property of solutions proved in [36] is based on showing 
a contraction property on the map $ acting on M =^A1i(C([— r, T], ]L^(r)) the space of laws 
of stochastic processes taking values in ]L^(r) ^=Ij\{T,E) defined, for any X G by: 

'Coir) + Jo ds(G{r, s,X,{r)) + A(r') dr' d7j{r, r', u)lEz[b{r, r',X,{r), 

+ A(r') dr' d^i{r, r', u) Mz[l3{r, r', X,{r), Z,+„(r'))] • dB,{r, r') 
+ J^dWs{r)g{r,s,X,{r)) , t>0 



C?(r) 



te[-T,0] 



(Zt) ={Xt) e M independent of (Xt), {Wt{-)) and •)) 



It can be shown that the unique solution of the mean-field equations is the limit of the recur- 
sion Xj^~^^{r) = ^{X^{r)) starting from any initial process X^{r). Let X^{r) a stochastic 

process whose law does not depend on r: X^{r)=X^{ro) for any (r,ro) G T. We can then 
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show that <l>(X^)^(r) do not depend on r either, since we have: 

<f{X%{r) = X°{r) + ds(G{r, s, X°{r)) + ^ A(r') dr' j' d^{r, r', u)E^o [5(r, r', X^r), ^°+„(r'))]^ 

+ / \{r')dr' r d^i{r,r',u) fw^z^Wy ,X'i{r),Z%^{r'))]- dB,{r,r') + f dW,{r)g{r,s,X°{r)) 
Jr J-T Jo Jo 

^X^{ro) + jys(G{ro,s,X°{ro)) + J^X{r')dr' J° dr^ir,/ ,u)]Ezo[b{r,r' , X°{ro), Z°+^{ro))]) 

\{r')dr' r dfi{r,r',u) f ]Ezo[/3{r,r' , X°{ro), Z°+^{ro))] ■ dB,{ro,r') + f dWs{ro)g{ro,s,X°{ro)) 

J-T Jo Jo 

= X°{ro) + ds(G{ro, s, X°{ro)) + E^o [B{r, X°{ro), 4)(ro))]) 

+ Mzo[H{r,X'^{ro),ZlJro))] + f dW,{ro)g{ro, s, X°^{ro)) 

Jo 

^X°{ro) + ds(G{ro,s,X°{ro))+^z[B{ro,X°{ro),Zl^{ro))]) 

+ Mzo[H{ro,X°,{ro),Z°{ro))] + f dW,{ro)g{ro,s,X°,{ro)) 

Jo 

= $(XO)t(ro) 

because of the assumption (6). The sequence of processes XJ^(r) is hence a sequence of 
processes whose law is independent of r by an immediate recursion, and so is the hmit. 
Hence the unique solution is spatially homogeneous in law. Moreover, a fixed point of 
the solution of the mean-field equation satisfies equation (7). 

The existence and uniqueness of solutions of equations (7) is more complicated to prove 
and necessitates different stochastic calculus arguments. It is outlined in Appendix A. □ 

A similar sufficient condition can be derived for polychronization in law is more intricate 
notationally. Sufficient conditions of the type of that of proposition 3 would involve defining 
a partition of T into different clusters and ensuring that the input received by any neuron in a 
particular cluster from the other clusters only depends on the cluster the neuron belongs to. 
Conceptually, this does not differ from the full synchronization in law condition, and would 
correspond to a new partition of the neural field F into subpopulations in which neurons 
are interchangeable. This condition is more delicate to write. Pathwise synchronization of 
mean-field neural fields, i.e. the actual convergence of the trajectories of the solution of the 
network equations for two different neurons is more complicated to ensure in our mean-field 
setting. Such a property would for instance involve a stochastic contraction property of 
the mean-field equations that would indeed ensure that the distance between two solutions 
starting from different initial conditions decreases towards zero. Such properties are not 
addressed here. 

The present section summarized a few theoretical results and gave a formal condition 
for spatially homogeneous solutions to exist. This very formal approach is now used to 
investigate a simple neural field model that is amenable to analytic treatment, firing-rates 
neuron models. We start by theoretically deriving the reduced equations in section 3 before 
turning our attention to their dynamics in section 4. 
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3. Mean-field equations for firing-rate models 

We revisit, from our probabilist point of view, the seminal works of Wilson and Cowan 
for finite-populations networks [20] and neuronal tissue [21], and the neural fields with prop- 
agation delays [46]. In the work of Wilson and Cowan, the authors emphasized the impor- 
tance of the presence of delays in the dynamics, a feature which has long been overlooked 
in the neuroscience community, but that has recently proved essential to shape cortical 
activity [5, 7, 47]. In that model, the state of each neuron is described by a scalar quan- 
tity representing the voltage of each neuron, assumed to have a linear intrinsic dynamics 
Ga{t^x) = —x/6a -\- Ia{t) (respectively G{r^t,x) = —x/0{r) -\- I{r,t) in the spatially ex- 
tended case), that interact through their mean firing rate which is given by a sigmoidal 
transform of the voltage variable. The delayed interactions are written as a sum over all 
neurons of a synaptic coefficient only depending on the populations the interacting neu- 
rons belong to and a sigmoidal transform of the pre-synaptic neuron (the one sending a 
current), baj{x^y) = Ja-fS^{y) (resp. b{r^r\x^y) = J{r^r')S{r' ^y)). The functions S^{x) 
and S{r,x) are assumed uniformly bounded and uniformly Lipschitz-continuous with re- 
spect to X. We enrich the model by taking into account intrinsic noise, i.e. that each 
neurons integrates a stochastic current driven by a Brownian motion, with a diffusion coef- 
ficient ga{t^x) = Xa{t) (resp. g{r^t) = A(r, t)). We also consider that the interconnection 
weights are noisy, as done in section 2, and specify the function /3c^^(x, y) = (Ja-fS^{y) (resp. 
I3{r^r' ^x^y) = (j{r^r')S{r' ^y)). These models clearly satisfy the assumptions of theorems 1 
and 2, and we use these to derive a stochastic version of their equations as the number of 
neurons tends to infinity. In order to simplify further our analysis, we will assume that the 
delay measures r]aj{u) and fia,-f{u) (respectively r]{r^r'^u) and fiir^r' ^u)) are Dirac mea- 
sures at fixed times t^^ (respectively at r{r^r')^ which generally can be considered to be 
Ik — ^'||r/c where c would the transport velocity in the axons). 

The equation of the dynamics of neuron i of population a in the network with N neurons 
reads, in the case of a finite number of populations, 

/ 1 P X 

^F^'^(t)= -^l^^'^(t)+/.(t) + ^J«,-^^5,(F^-'^(t-r.,))Ut+ 

\ 7=1 ^ j = l / 

Xa{t)dWi + ^ a^, J2 S.iV^^'^it - ra,))j dB^ (8) 
and in the case of spatially extended networks: 



= (^-^^'•'^W+^(''-*)+E J^^-^^f) j^J2S{r„V^^''{t-r{ra,r,)))j dt+ 

P{N) / 1 \ 

A{r^,t)dWi+ air^,r^)\^—J2sir^,V^'''it-Tir^,r^)))j dBl{r 



K). (9) 



Note that we will also consider neural fields composed of different layers (for instance exci- 
tatory and inhibitory neurons). The above theoretical analysis will readily extend to that 
case. 
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In these models, Oa (resp. 0{ra)) are the characteristic time of return to rest of neurons 
in population a if it is not submitted to any input, and these times are assumed to be 
lowerbounded by a strictly positive value Om- The current Ia{t) (resp. I{ra,t)) corresponds 
to a deterministic input current fed to the neuron, J^^ (resp. Jivo^^r^)) is the synaptic 
weight of the incoming connections from neurons of population 7 to neurons of population 
a and Gcx^ (resp. <j(ra,r^)) its standard deviation, \a{t) (resp. A(ra,t)) the standard 
deviation of the external noise, and the functions S^{-) (resp. 5'(r^, •)) are smooth sigmoidal 
transforms, assumed bounded with bounded derivative. 

3.1. Finite-population Wilson and Cowan Equations 
As an application of theorem 2, we have: 

Proposition 4. Let T > a fixed time. Assume that the initial condition on the network 
are chaotic, i.e. independent and identically distributed for neurons in the same population. 
The the state V^'^ of neuron in population a, solution of equation (8), converges in law as 
N goes to infinity towards the process Va unique solution of the mean-field equation: 



dVa{t) 



-^v^{t) + i^{t) + Ja^ns^w - 



^=1 



dt 



+ K(t)dW^ + ^ CJo^^nS^W - Ta^WBf (10) 

where the processes (W^t) and B^^ for a,7 G {l...P}^ are independent standard real 
Brownian motions, and the propagation of chaos property applies. 

Moreover, if = {V^)a=i---p ^ M'^{Cr) is a P-dimensional Gaussian process with in- 
dependent coordinates, the solution of the mean- field equations (10) with initial conditions 

are Gaussian processes for all time, with mean denoted fi{t) = {/J^a{t))a=i...p (^nd vari- 
ance v{t) = {va{t))a=i...p- Let fi3{x^y) denote the expectation of Sp{U) for U a Gaussian 
random variable of mean x and variance y. We have: 



f^a{t) = --X-I^a{t) ^y^^ Ja(3f(3(l^(3{t - ra(3),V(3{t - Tafd)) +^a(^) a= 1...P 

2 ^ 

Va{t) = -y Va,(t) ^Y(jl^fl{^^(t - To,p),Vp(t - To,^)) +A^(t) a= 1...P 



(11) 



with initial condition ii^it) = lEj[V^{t)] and v^it) = ]Ej[{V^{t) - /ia(t))^] for t e [-r,0]. In 
equation (11); the dot denotes the differential with respect to time. 

Proof. It is indeed easy to show that the network equations are of the form of equation (1) 
with smooth parameters that satisfy the assumptions (HI)- (H4). The first part of the 
proposition is a direct application of theorem 1. The P equations (10), which are P implicit 
stochastic differential equations, describe the asymptotic behavior of the network. The 
unique solution of the mean- field equations (10) with initial condition can be written in 



Stochastic Neural Fields Dynamics 



14 



an implicit form as: 

+ V f et;a^^nS^{V^{s-r^^mdBf ^ f et; \^{s) dW^) . (12) 

13=1 

It is clear from this formulation that the righthand side is a Gaussian process as the sum of 
a deterministic function and of a sum of stochastic integrals of deterministic functions with 
respect to Brownian motions, and hence so is Va{t)^ - Its distribution is hence characterized 
by its mean and covariance functions. The term ^[Sj3{Vj3{s — Tc^/?))], because of the Gaussian 
nature of Vg, only depends on lJi^{s — Ta/s) and Vi3{s — Ta/s)^ and is denoted by f(3{ii^{s — 
'Ta^)i Vi3{s — Ta/s))' Taking the expectation of both sides of the equality (12), we obtain the 
equation satisfied by the mean of the process jj^ait) = El[Va(^)]- 



fia{t) = e ^« |/ia(0) + J //3(/i/3(s-rc,/3),V/3(s-rc,/5)) +/«(<§) I (is I . 

Taking the variance of both sides of the equality (12), we obtain the following equation: 

Va{t) = e-^ |^^c.(0) + ^ |^^^^'/3/|(M^(5-^c.^),^/5(5-rc.^)) + A^(s)j ds^ . 

These two formulae are exactly an integral representation of equations (11) which concludes 
the proof. □ 

Remark. 

• The variance by itself does not characterize the law of the process V, it is necessary to describe 
the covariance matrix function Co,p{ti,t2) for ti and t2 in R^*. It is easily shown that for 
a ^ (3, this correlation is null, and moreover: 

ti+t2 ^ /-tiAts 2. ^ 

Ccc.(ti,t2) = e ^« v'^{0)^ e^c.\^{sYds 

Jo 

for ti,t2 G R^*; hence only depends on the parameters of the system and on the functions 
lip{t) andvf3{t). The description of the solution given by equations (11) is hence sufficient to 
fully characterize the solution of the mean-field equations (10). 



^This property can also be proved by using the classical characterization of the solution of the mean-field 
equation as the limit of the iteration of the map <l> as proved in [36]. In that case, similarly to the proof 
provided in [39] for non-delayed finite-populations equations, V is defined the limit of a sequence of Gaussian 
processes, hence Gaussian itself. 
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• In the case where the sigmoidal transforms are of the form Sa{x) = erf{gaX + ha), a simple 
calculation based on a change of variable shows that the functions Va) involved in the 

reduced mean-field equations (11) take the simple form (see [39] for the calculations): 



This simple expression motivates the choice of erf sigmoidal functions for further develop- 



We can hence reduce the complex mean-field equation (10) into the simpler system 
of coupled delayed differential equations (DDEs) (11), which allows a simple analysis of 
the solutions and the dependency of the behavior upon the parameters. Even if we know 
that there exists a unique solution to the mean-field equations, and that necessarily the 
solution starting from Gaussian chaotic initial condition is Gaussian with mean and standard 
deviation satisfying equations (11), it is necessary to show that these latter equations are 
well-posed and that they have a unique solution. Let us denote by C the Banach space of 
continuous functions mapping [— r, 0] into E'^^ endowed with the topology of the uniform 
convergence. Following Hale and Lunel [48], we consider the moment equations (11) as 
ordinary differential equations on C. We have the following: 

Theorem 5. Under the assumptions that: 

• t^ Ia{t) and t ^ Xa{t) are continuous, 

• Xa{t) > Ao > for any t > and a G {1, • • • , P}, 

• the sigmoid functions have bounded derivatives 

• The initial conditions on the variance v^, e C([-r,0],R) satisfies vl{t) > vo > for 
t e [-r,0] and a G {1, • • • ,P}. 

there exists a unique solution to the moment equations (11) starting from an initial condition 
111 e C([-r,0],R) andv^. 

Proof. This theorem is a simple application of theorems [48, Thm 2.1 and 2.3]. These 
results ensure existence and uniqueness of solutions as soon as the vector field is Lipschitz- 
continuous in C and continuous in time. Thanks to uniform lowerbound of the functions 
and of the initial condition on the variances v^, we clearly have v{t) > Vm = m.m{vo^ X^Om/'^) 
(we recall that Ojn is the strictly positive lower bound of the characteristic times ^q,). We 
have: 



Thanks to the properties of the integral term and in particular the fast convergence towards 
of the exponential term ensuring existence of moments of any order of the Gaussian, it is 
straightforward to show that derivative with respect to x and y read: 




ments. 




dy 



dx 




dz 



dz 



(13) 
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which are both uniformly bounded on R x [vm^ oo) thanks to the fact that the sigmoids have 
bounded derivatives and that y > vq^ since we have: 





df^ix,y) 




dx 


j 


df^{x,y) 




dy 



ll'^illc 



< 



2^ jR 



/r 



dz 



\\s'\ 



Jo 



^/2 



dz 



WS'J 



This property ensures global Lipschitz continuity of the vector field. The differentiability of 
the vector field with respect to time readily follows from the differentiability of the input and 
standard deviation functions. We are hence is the application domain of theorems [48, Thm 
2.1 and 2.3] and we have existence and uniqueness of solutions for the reduced mean-field 
equations (11). □ 



Remark. Note that the conditions of theorem 5 are not optimal. 
Sa{x) = erf{gaX + ha), we already remarked that we had 



For instance in the case 



U{fJ^,v) erf 



and this function is globally Lipschitz continuous on R x R^^ as one can easily check by computing 
exactly the derivatives of this function with respect to /j and v. The condition on the uniform 
lowerbound of Xa (t) in theorem 5 can hence be omitted when one chooses erf sigmoidal transforms 
for instance. 



3.2. Spatially extended neural fields 

We now consider the infinite-population case of firing-rate neuron models. We assume 
that we are a priori given a connectivity map with mean J(r, r') and standard deviation 
cF{r.,r')^ a density A(r) such that the neural populations are drawn independently according 
to this law (after renormalization) as the number of neurons increases. Using theorem 2, we 
can state the following proposition: 

Proposition 6. Let T > a fixed time. The process T/*'^ for i in population a at location 
ra ^ solution of equation (9) with chaotic initial conditions, converges in law towards the 
process V{ra) where {Vt{^)) ^-^ the unique solution of the mean-field equation: 



dVtira) 



^Vt{ra) + I{ra,t) + f \{r')dr'J{ra,r')nS{r',Vt.r(r^y){r'))] 
^a) Jr 



0{ra) 



dt 



+ K{r^,t)dWt{ro) + J^a{r^,r'MS{r',Vt_r(r^y){r'))]dBt{r^,r')X{r')dr' (14) 



where the processes {Wt{r)) and {Bt{r^r')) are independent cylindrical Brownian motions, 
and the propagation of chaos applies. 

Moreover, if the initial condition F^(r) G A^^(C([— r, 0]), ]L^(r)) is a cylindrical Gaus- 
sian process, the solution of the mean- field equations (10) with initial conditions V'^(r) is 
Gaussian for all time. Denoting by /J^{r^t) its mean, by v{r,t) its variance, and by f{r,x,y) 
the expectation of S{r, U) for U a Gaussian random variable of mean x and variance y. We 
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have: 
dt 



(r, t) = - ^/i(r, t) + / X{r')dr'J{r, r')f{r, fi{r\ t - r(r, /)), v{r\ t - r(r, /))) + /(r, 



^(r, = Kr-, + Hr'fdr'a{r, r'ff{r, fi{r', t - T{r, r')), v{r' , t - T{r, r'))f + K\r, t) 

(15) 

with initial condition ii{r,t) = E [l^t°(r)] and v{r,t) = E[(y°(r,t) - /i(r,t))^] /or t G [-r, 0] 
anti r G r. 

Proof. This proof is essentially identical to that of proposition 4. Indeed, it again appears 
clearly that the network equations are of the form of equation (4) with parameters satisfying 
the assumption of section 2.2, and hence that theorem 2 applies. The unique solution of the 
mean-field equations (10) with initial condition satisfies the implicit equation: 

(16) 



and is hence necessarily Gaussian since the righthand side obviously is. Its mean, denoted 
by /i(r, t), and variance v{r^t)^ satisfy the equations: 



/i(r,t) =e ^(o(^^(r,0) + /oe^^)(/(r,5) 

+ \{r')dr' J{r^ r')f{r\ /J.{r\ s — T{r^ r')), v{r' ^ s — T{r^ r')))^ds^ 
v{r,t) =e"^^v(r,0) + (^A(r, s)^ 

+ \\r')dr'cj\r, r')p{r', /i(r^ s - r{r, r')), v{r\ s - r{r, n)))ds^ , 

which concludes the proof. □ 

Remark. Here again, mean and variance are enough to define the law of the stochastic process, 
since the covariance matrix between V{r,ti) and V{r\t2), C{r,r ,ti,t2) for ti and t2 in R^* and 
{r,r) G r only depends on these two quantities. Indeed, for r ^ r' , the covariance is null because 
of the cylindrical nature of the Brownian motions involved and of the initial conditions. It is easily 
shown that for r — r' : 

rtiAt2 



C(r,r,ti,t2) = e"^^(^^^;"(0)+ / e^A{r,sfds 

Jo 

/*i/\t2 2s r 
e^M / X(r)^dra{r,r)^f^{r,fi{r,s — r{r,r)),v{r,s — r{r,r))i 



The description of the solution given by equations (15) is hence sufficient to fully characterize the 
solution of the mean-field equations (14). 
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We therefore describe in this formahsm the stochastic dynamics of a neural field through 
two coupled integro-differential equations, instead of just one in the classical cases: the 
standard deviation of the process is now incorporated to the dynamics of the mean activity, 
and in the non-noisy case, correspond exactly to more classical neural field equations [46, 



Similarly to what was done in the finite-population case, it remains to show that re- 
duced system on the moments (14) is well-posed, to ensure that these equations uniquely 
characterize the mean-field equations. This is the case under the following conditions: 

Theorem 7. For the sake of simplicity, we assume here that the density function A(r) is 
upperbounded by a constant A. Under the assumptions that: 

• J is square integrable with respect to Lebesgue's measure on , i.e. belongs to ]L^(r^, R)^ 

• cr^ is square integrable with respect to Lebesgue's measure on , 

• the external current I{r^ t) is a bounded, continuous functions of time taking values in 



• the external noise K^{r^t) is a continuous functions of time taking values in ]L^(r,R)^ 
and is uniformly lowerbounded by a strictly positive constant: A{r,t)'^ > Aq > for all 



• The derivative of 5'(r, x) with respect to its second variable is uniformly bounded, 

then for any initial condition ii{r,t) G C([-r, 0], ^^(r, R)) and v{r,t) G C([-r, 0], ^^(r)) 
uniformly lowerbounded by a quantity vq > 0, there exists a unique solution to the moments 
mean-field equations (14) which moreover belongs to C([— r, T], L^(r, R^)). 

Proof. The moment mean-field equations (14) constitute a dynamical system in the Banach 
spaces of functions of F with values in R^. It is well known that the space of functions in 



is a Banach space. We will show the existence and uniqueness of solutions in this space. We 
further define the norm up to time t > of two elements of B by: 



Let us start by ensuring that any possible solution is bounded in this space. We recall 
that: 



49, 50]. 



]L2(r,R) 



(r,t) G r X R+, 



S = C([-r,T],L2(r,R2)) endowed with the norm: 
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and hence we have: 



+ [ dr\ f ds f X{r')dr'J{r,r')f{r',ij{r',s-T{r,r')),v{r',s-T{r,r')))\' 
Jr Jo Jr 



< 3 ||E [V°{-)] \\l^r,K)+T^ sup \\I{;s)\\l,^r,K) 

\ *e[o,T] 

+ rA(r) ^ dr ds ^ \{r')dr'\J{r, r')|'||/||L) 

< 3 (||E [T/o°(-)] \\iHr,n)+T^ sup ||/(-, s)||^(r.R) + T^A A(r)|| J(r, r')||i2(r2,R)| 

\ te[o,T] 

where ||/||oo is the uniform upperbound of /(r, /i,!;) in R, which is smaher or equal to the 
uniform supremum of the function x S{r,x) (which exists by assumption). The same 
types of calculations allow proving that: 



\H-MlHr,n)<^[\H-M\lHr,n)+T' ^sn^^^^^ 

These bounds do not depend upon time t and hence prove that any solution of the 
moment mean-field equations have bounded norms in the space B. 

Routine methods for this type of infinite-dimensional systems ensure existence and 
uniqueness of solutions as soon as the vector field of the equation is Lipschitz-continuous for 
this norm. In our case, similarly to what was done in the proof of theorem 5, it is easy to 
show under the assumptions of the theorem that the standard deviation v{r^t) is uniformly 
lowerbounded by a minimal strictly positive value v^n = min('Uo, and hence ensur- 

ing the global uniform in r Lipschitz-continuity of the function {ji^v) ^ f{r^ii^v). Let us 
define for {(i^v) G B the transformation ^{ji^v) taking values in B and defined by: 



v 



I oo 



^(O r/i(0, r) + /q e ^(-) (^/(r, s) + \{r')dr' J{r, r')f{r', ii{r' , s - T{r, r')),v{r', s - T{r, r')))^ ds^ 
^ ^^(0, r) + /q (A^(r, s) + X{r')'^dr'(j'^{r, r')p{r', /i{r\ s - r{r, r')), v{r\ s - T{r, r')))^ds^ 



It is clear that any solution of the moment equations are fixed points of ^ and reciprocally. 
Since (S, || • II23) is a Banach space, showing existence and uniqueness of solutions, i.e. of 
fixed points of <l>, amounts showing a contraction property on ^. First of all, similarly 
to what was done to show that any solutions of the moment equations were bounded in 
S, it is very easy to show that for any {fi^v) G S, we have ^{ji^v) G B. We use the 
classical iteration method to show existence and uniqueness of fixed point. To this end, we 
fix ip^ = G B arbitrarily and define the sequence cp^ = (/i"', '^"')nGN iteratively by 

setting 1'^+^) = ^{jii'^ ^v'^). We recall that / is Lipschitz-continuous as shown in the 

proof of theorem 5, and we denote by L the uniform Lipschitz constant of fir^x^y) in its 
two last variables. The function p{r^x^y) is hence also uniformly Lipschitz-continuous in 
its two last variables with the Lipschitz constant 2||/||oo^- Let us now show that the vector 
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field $ is Lipschitz-continuous on B. Let us fix cpi = (/ii, t'l) and (p2 = (/i2 5 '^2) two elements 
of S. We have: 

Dt{(pi,(p2)= sup \ dr X{r')dr'J{r,r')(f{r\fii{r',s-T{r,r')),vi{r\s-T{r,r'))) 
se[-T,t] I Jr ^0 ^ 

X 2 

- /(/, /i2(r', s - T(r, /)), ^2(r', s - T(r, /))) jc^s 
+ / / / >^\r')dr'a\r, /) . - r(r, . - r(r, /))) 

- /^(r', /i2(r', s - T{r, r')),V2{r' , s - T(r, r')))^ds | 

The two terms of the righthand side are treated similarly, let us hence deal with the first 
one. We have: 

J J J ^{'^')dr J{r^r)(^f{r\ jj.i{r\ s — r{r^r))^vi{r\ s — r{r,r))) 

- /(r', /i2(r', s - T{r, r')),V2{r', s - r{r, r')))jds 

<T f dr f\j \{r')dr'J{ry)(^f{r',^i^{r',s-T{ry))M 
«/ r «^ «>' r 

- f{r', ii2{r', s - T{r, r')),V2{r', s - T{r, r')))jds 

- /(^', M2(r', s - T{r, r')),V2{r\ s - T{r, r')))^ j ds 
<2L^T J dr (^j \(r')dr'J{r,r'f^(^j \{r')dr'\^ix{r' , s - T{r,r')) - ii2{r' , s - T{r,r'))\^ 

+ \vi{r' ,s - T{r,r')) - V2{r' ,s - T{r,r')))\^^ds 

<2A^L^T j^dr j'^ dr'J{r, ||^i ^ - r(r, /)) -^'.s- r{r, /)) ||^(r,R^) 

<2A'L' TWmrxr) f ^.(^1, ^2) ds 

Jo 

Similarly, the second term is upperbounded by 8 ||/||^ ||cr^ ll^s^rxr) Jo -^s {^1:^2) ds. 
These two bounds do not depend upon time, hence we have: 

^2) <K' Ds{(pi,(p2) ds 
Jo 

with iC' = 2 ( 1 1 J I (r X r) + 4 I 1 1 i2 (r X r) ) • Routine methods allow showing that : 
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and hence ||(/:'^ — ^^~"^||i3 < \/ {K^t)'^ /n\\\(p^ — (P^\\b^ readily implying that cp^ is a Cauchy 
sequence in the complete space B and hence converging in B towards a fixed point of 
Uniqueness of the solution is also classically deduced from the above inequality. □ 

The well-posedness of the moment equations allows studying the solutions of the stochas- 
tic mean-field equations through the analysis of the dynamical system (15). 

4. Analysis of the solutions 

Now that mean-field equations for firing-rate models were derived, we numerically and 
analytically analyze their solutions. We investigate separately the effect of delays on the 
solutions in a finite-population network, and spatially extended neural fields equations with 
or without delays. The specific effects of propagation delays in neural fields is not analyzed 
here. 

4-1- Finite-populations networks 

We analyze in this section the effect of delays in the dynamics of finite-populations 
neural fields. To fix ideas, we start by a numerical analysis of a simple 2-populations model 
composed of one excitatory and one inhibitory population, with synaptic weights: 

^ / 15 -12 \ 

and constant delays r. The time constants Oi and O2 are set to 1, and the noise on the 
synaptic connectivity ctq,^ are all null. The currents are fixed to Ii = 0, /2 = —3, and both 
populations have a common noise intensity A. The bifurcation diagram of the system (11) 
as a function of the noise intensity A and the delay r is produced in Figure 2. Two branches 
of fixed points are identified, and as delays are increased, one of the branches of fixed 
points undergoes a cascade of Hopf bifurcations. The different curves of Hopf bifurcations 
accumulate on the first curve as delays are increased, on a parameter region where the related 
fixed point is unstable (see Figure 2(e)). These Hopf are not related to eigenvalues having the 
largest real part, and hence has no effect on the number or stability of fixed point. These Hopf 
bifurcations are associated with unstable limit cycles, the accumulation of which produce a 
complex landscape resulting in very irregular transient behaviors. When fixing a value for 
the delays to r = 0.5 (green line of figure 2(a)), we observe in this diagram different ranges of 
parameter values corresponding to different asymptotic behaviors: stationary solution (blue 
region), bistability between a stationary and a periodic solution (yellow region), and periodic 
solutions (orange region, see Fig. 2(e)). The oscillatory region corresponds to a precise 
synchronization of all neurons (see figure 2(f)), a highly relevant phenomenon in neural 
systems sometimes related to pathologies such as epilepsy. Fixing r = 5 corresponds to a 
case where the system displays seven Hopf bifurcations. The codimension two bifurcation 
diagram of the system as a function of the input current Ii and the noise level A is given 
in Figure 2(c). We observe that the saddle- node bifurcation forms a cusp, and on one 
of the branch of the codimension two saddle-node bifurcation curve appears a Bogdanov- 
Takens and two degenerate Bodganov-Takens bifurcations. These degenerate bifurcations 
correspond to the tangential pairing of two Hopf bifurcations merging with the saddle- 
node manifold. These are non-generic bifurcation but appear structurally stable in these 
equations. The precise study of the unfolding of these bifurcations is not in the scope 
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of the present paper, but the unfolding of that bifurcations would be helpful to further 
understand the local behavior of the system around these singular points. An interesting 
feature of the system is that delays create complex phenomena only in a restricted region 
of parameters corresponding to small noise values: for noise large enough, no instability 
is found, illustrating quantitatively the stabilization effect of the noise. It also illustrates 
the strong impact of noise in shaping the qualitative form of the solution of the mean-field 
equations: for small noise, a stationary behavior is found, for intermediate values a periodic 
orbit is found corresponding to the synchronization of different neurons in the network (see 
Figure 2(f)) and for larger values of noise, a stationary behavior is again obtained. 

The previous example was particularly interesting because noise induced the presence of 
regular oscillations. However, its particular dynamics and parameters do not allow further 
mathematical investigation because the characterization of fixed point was not analytically 
explicit and hence it was hard to characterize the local bifurcations. In order to analytically 
investigate this kind of system and in particular account for the cascade of Hopf bifurcations 
numerically exhibited in the previous example, we now simplify the problem by considering 
a similar neural network model with connectivity matrix: 

and input currents /i = and I2 = —I. In this case, /ii = /i2 = 0, = = A^/2 is a 
fixed point of the equation. The characteristic matrix governing the linear stability of this 
solution reads: 

whose eigenvalues (the caracteristic roots) are: 

v± = -(C + 1) + ^ =e-<"(l ± i) 

with = —1, and the characteristic equation, defined as A(^) = det{A{C,)), vanishes for 
values of C such that at least one of the characteristic roots vanishes, i.e. for C such that 

-(C+l) + ^ ^ =e-<"(l±i)=0 (18) 

v/27r(l+ff2A2/2) ' ^ ' 

This equation can be solved using the complex branches of Lambert's {Wk)ke'E functions 
(see e.g. [53]). Simple algebra yields the following formula for the characteristic roots of the 
system: 

Ci = -1 + -Wk I ^ =re^l ± i) 1 . (19) 

^ r ' VV2^(1+^'AV2) ^ V ^ ^ 

The stability of the solution is governed by the sign of the real part of the uppermost 
eigenvalue. There can be an infinite number of characteristic roots of the equations, but for 
any real there exists a finite number such that 5R(C) > 1^0 - The eigenvalue with largest 
real part is given by the real branch Wq, and if the argument has a real part greater than 
—e~^ the root is unique. If this is not the case, two eigenvalues have the same real part, 
namely the ones corresponding to/c = Oor/c = — 1. 
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T 10 




(a) Codimension 2 bifurcations in (b) 3 first Hopf bifurcations, (c) Codimension 2 bifurcations 
(A,r) large delays in (/i, A), r = 5 




(d) solution, t l~J^i(t) for r 
10, A = 0.1 





(e) Codimension 1 bifurcation (f) 100 trajectories for the net- 
diagram, r = 0.5 work equations, A = 1.5, r = 0.5 



Figure 2: Dynamics and bifurcations of the mean-field equations (11). (a): Codimension 2 bifurcation 
diagram as a function of the noise intensity A and the delay r: saddle-node bifurcations (blue line) and 
cascade of Hopf bifurcation (pink curves) that all have a common vertical asymptote (b). (c): Codimension 
2 bifurcation diagram as a function of Ii and A for r = 0.5 (green line in diagram (a)): 2 degenerate 
bifurcations appear, corresponding to the tangential merging of two Hopf with the saddle-node bifurcation, 
and one Bogdanov- Taken bifurcation, (d) Codimension 1 bifurcation diagram as a function of the noise 
intensity A and Ii = 0, segmented into three zones: blue for stationary solution, yellow for bistability 
between a stationary/ periodic solutions, and orange: periodic behavior. In the periodic region, all neurons 
in the network are synchronized (2(f): blue (resp. red): 50 trajectories from neurons of the excitatory 
(resp. inhibitory) population 1 (resp. 2)). The diagrams were obtained using DDE-BIFTOOL [51, 52] and 
a specific code for the network equations. 
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We observe that the argument of the Lambert function in the expression (19) has a 
modulus that decreases towards as A or r are increased towards infinity. For fixed values 
of r, the rightmost eigenvalue is given by /c = 0, we hence observe that as A is increased, the 
value of Wo decreases towards 0, and there exists a value A(r) such that for any A > A*(r) the 
fixed point is stable: noise has a stabilizing effect on this fixed point. This stabilization 
appears through a Hopf bifurcation, and periodic behaviors are found for A < A*(r), as 
shown in the bifurcation diagram 3(e) for a fixed value of the delays, r = 0.5. For fixed 
values of A, as r is increased, the real parts of the eigenvalues increase and might switch 
from positive to negative (see figure 3(a)). As delays are increased, a growing number of 
eigenvalues have positive real part, which accounts for the more complex phenomena of 
creation and destruction of limit cycles and for the increasing number of fixed points as a 
function of A and r. 

Let us now investigate more thoroughly the location of possible Hopf bifurcations. These 
are necessarily points where the real part of the characteristic roots vanishes. Following [48, 
54], we use the fact that a necessary condition for the existence of a Hopf bifurcation is 
the existence of a purely imaginary characteristic root ( = iuj. This condition, plugged into 
formula (18) yields the condition: 

-(iuj + 1) = ^ =e-^^^(l ± i) (20) 

which, taking the squared modulus of these imaginary numbers, give the equality: 

1 -\- uj — 



-7r + ff^(l-7rAV2) 



This quantity is necessary positive, hence Hopf bifurcations only arise for values of the noise 
intensity satisfying the inequality: 

This inequality precisely corresponds to the vertical asymptote: for any value of the noise 
intensity greater than A*, no Hopf bifurcation is possible in the system. This property 
indicates that necessarily, for noise intensities greater than A*, the fixed point (0,A^/2) is 
stable. Indeed, we already mentioned that for large noise values that fixed point is stable. It 
can gain stability only if the real part of one of the characteristic roots crosses the imaginary 
axis, corresponding either to purely imaginary roots or null roots, and the previous inequality 
showed that such roots do not exist for parameters A > A*. 

Under the condition A < A* , we can perfectly characterized Hopf bifurcations by equating 
the argument of both sides of equality (20). We obtain that way: 

— arctan(a;) ± | + 2 /ctt 



for /c G This relationship can be written in closed form as a function of the parameters 
using the expression of u obtained in equation (21). The different curves of Hopf bifurcations 
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Figure 3: Analysis of the characteristic roots of the delayed neural field equations related to the equilibrium 
point fii = IJ'2 = and vi = V2 = Results are similar to those of Fig. 2: as delays are increased, 

several Hopf bifurcations arise and accumulate around the same value, (a): Shape of the Lambert function 
X ^{Wo(x(l-\-i))), (b): cascade of Hopf bifurcations and (c): locus of the Hopf bifurcations for the 15 first 
characteristic root continuated for very large values of r (d): Transient regime for r = 5 and (e): Bifurcation 
diagram as a function of A for r = 0.5. 



are plotted in Figure 3. We observe the exact same phenomenon of cascade, accumulation 
and pairing of Hopf bifurcations as delays are increased. The accumulation is related to the 
fact that the value of A related to a Hopf bifurcation increases as a function of r, and that 
the value of A is upperbounded. This case again is characterized, for large values of r, by the 
presence of irregular transient behaviors corresponding to the very complex landscape of the 
phase plane, as displayed in Fig. 3(d). We chose for instance to display the transient solution 
for r = 4, a case where the delays are small enough so that we can resolve the presence 
of different limit cycles trapping the solution transiently. This case hence makes explicit 
the dependence on noise levels of qualitative behaviors of the system for finite populations 
networks. We now extend the study to infinite populations networks in a neural field setting. 

4.2. Neural Fields Dynamics 

We now turn our attention to the study of the dynamics of the neural field equations 
derived in section 3.2. We start by an analytic study of the dynamics of a single layer 
neural field with delays with a particular focus on Turing-Hopf bifurcations from spatially 
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homogeneous states, before addressing numerically the same questions in the more complex 
and more relevant dynamics of a two-layers neural field composed of an excitatory and an 
inhibitory layers. 

4.2.1. A single layer neural field: Turing- Hopf instabilities and noise 

In this section, we consider the case of a single-layer neural field, whose mean and 
standard deviation satisfy equations (15). In order to analyze the dynamics of the solutions, 
we chose a particular set of parameters that allow analytical studies. We consider a one- 
dimensional periodic neural field on F = distributed homogeneously (i.e. X{x) = 1), with 
connectivity functions J{r^r') and a{r^r') only depending on |r — r'|, 0{r) constant, and 
denote with a slight abuse of notation the constant value of the function ^(r), J{r^r') = 
J{r — r') (resp. cr(r, r') = cr(r — r')). We denote by J (resp. a^) the integral J{r.,r')dr' 
(resp. cr^(r, r')dr') assumed finite (these obviously do not depend on r). Taking 5'(r, x) = 
erf(^x), the function f{r^x^y) is equal, as already stated, to eTf{gx / y^^Y^^^^) . Hence 

^0 =V(^7 ^) does not depend on v. Let us further set / = —JFq. In that case, /i(r, t) = 
for any (r, t) G F x R+ is a solution of the mean equation whatever the standard deviation 
V is, and this variable v is solution of the equation: 

g(r,i) = -^t;(r,t) + a2Fo' + A2(r,t). 

Let us further assume that A(r, t) is constant. Then the spatially homogeneous, constant in 
time function: 

is solution of the neural field equations. 

It is clear that since the domain is periodic and both the connectivity and the delay 
functions only depend on |r — r'|, the related neural field satisfies the full synchronization 
property of proposition 3. The spatially homogeneous state (or synchronized state) satisfies 
the equations: 

^ = - f + /o Jir)f{pi{t - r(r)), v{t - r(r))) dr + I 
t = -¥ + /o '^\r)f{t,{t - T{r)), v{t - T{r))) dr + 

which is a delayed differential equation that is formally equivalent to a differential equation 
with delays distributed on the interval [min^ r(r), max^ r(r)] with the density J(r). The 
spatially homogeneous equilibria (/i, v) are hence solution of the equations: 

f-f + + ^ =0 

\-^-i + ~a^f{ji,v)+K^ =0 

and are obviously the same equilibria as for the non-delayed spatially homogeneous equa- 
tions: 

The bifurcation diagram of this non-delayed system as a function of the slope of the noise 
A is produced in Figure 4. If the slope g is small enough (precisely g < Y^27r/j7), the fixed 
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point (/io^'^o) is the unique fixed point and is always stable whatever the noise parameters 
are. For larger values of g (in Fig. 4 we chose ^ = 3), the system presents two distinct regimes 
as a function of the noise parameters: for small values of the parameters a and A, the system 
presents three fixed points, (/io^'^o) that is unstable and two distinct fixed points denoted 
and (/i2,'^2) which are stable, and that hence constitute spatially homogeneous and 
constant in time solutions of the mean-field equations. The system undergoes a pitchfork 
bifurcation on the manifold defined by ct^Fq + = —l/g'^ + j7/27r, and on this line the 
two fixed points (/ii, '^i) and (/i2, '^2) collapse with (/io, '^0) and disappear. For a and A such 
that o-'^Fq + A^ > + J'/'^tt system is left with a unique equilibrium, {jhq^ vq)^ which is 

now stable. Hence again appears a stabilization effect of noise on the spatially homogeneous 
fixed point {jj^o^vq). 

Let us now return to the full spatially extended system with delays, and address the 
problem of pattern formation in these equations beyond a Turing instability at the spa- 
tially homogeneous steady state (/io,'^o)- To this purpose, we analyze the linear stability of 
this fixed point (see e.g. [55, 56, 46]), which amounts characterizing the eigenvalues of the 
linearized equations around this fixed point. Since f{x^y) = eTf{gx/y^l + g'^y)^ we clearly 
have: 

■ a/ I g 1 d^f T^f 

22 

and hence the linearized equations around {jj^o^vq) read: 



= 



OA 

at 

OB 

dt 



(r, t) = - ^ A(r, t) + J{r' - r)A(r^ t - t{t' - r)) dr' 

(r, t) = -lB{r, t) + 2Fo {r' - r) A(r^ t - r(r' - r)) dr' ' 



We recall that the integral operators are convolutions on S^, these can hence be diago- 
nalized on the Fourier basis. Let us consider perturbations of the equilibrium of the form 
Ay^k(r^t) = 5R(e^^+^^^^) with v = icj + and leave B(r^t) unspecified. Let us denote by 
{ak(h>)) and (hki^')) the Fourier coefficients of the functions J(r)e~^^^^^ and (j{r)e~^^^^^\ 

{ak{v) = J(rOe"^^^^'^e"^^^^^' dr' 
[bk{iy) = j^(J^{r')e-^<'^'^e-^''^^'^' dr' ' 

The functions Ajy^j^ are eigenfunctions for the first equation the linearized system provided 
that u satisfies the relationship: 

^ = ^/e = -^F^ak, 

defining a dispersion relationship. The characteristic roots of the linearized equations is 
hence composed of the eigenvalues of the matrix: 

2FoF(^6fc -2/0 



which are exactly {—2/0, —1/0 -\- F^ak-, k G N}. In particular, this shows that no instability 
can occur on the standard deviation equation, and the whole stability of the homogeneous 



Stochastic Neural Fields Dynamics 



28 



fixed point only depends on tiie Fourier coefficients of the deterministic connectivity function 
multiplied by the exponential of the delay function, J{r)e~^''^\ and on the quantity Fq which 
depends on the different parameters of the system. Explicitly, in the original parameters, the 
spectrum is hence composed of the eigenvalues {— 1, — ^ + a/e . ^ k G N} with vq = 

|(<jFq + A^). Let aM be the largest Fourier coefficient of J: aM = max/c^^ cik- By our inte- 
grability assumption on J, this maximum necessarily exists because of Parcheval-Plancherel 
theorem. The solution (0,'Uo) is hence linearly stable as soon as —1/0 + aMg/V^ < 0. 

This particular form of the spectrum of the linearized operator quantifies directly a 
stabilization effect of the noise already mentioned. Indeed, let us assume that 6 and g are 
fixed. If — 1/^ + ag /y/27r < 0, then all the eigenvalues are negative whatever vq and hence 
the solution (0, vq) is stable whatever the noise connectivity matrix a and the additive noise 
A. If now —1/0 + ag/\^27T > 0, then for A and a small, the fixed point (O^vq) is unstable. 
When cr or A are increase, the fixed point will gain stability, since the maximal eigenvalue 
tends to —1/6 when vq goes to infinity^. 

An instability occurs when the characteristic roots such that u has a positive real part. 
A Turing bifurcation point is defined by the fact that there exists an integer k such that 
= 0. It is said to be static if at this point ^{lyk) = 0, and dynamic if ^{i^k) = cj^ 7^ 0. 
In that latter case, the instability is called Turing-Hopf bifurcation, and generates a global 
pattern with wavenumber k moving coherently at speed Uk/k as a periodic wavetrain. If 
the maximum of Xk is reached for /c = 0, another homogeneous is excited. 

In order to investigate the presence of Turing-Hopf instabilities, we specify further our 
system. We choose an exponential connectivity function J(r) = e"'^'/^ for some 5 > 0, 
and a delay taking into account both propagation (or axonal) delays corresponding to the 
transport of the information along the axons at a constant finite speed c, and constant delays 
accounting for the transmission time of the spike through the synaptic machinery. In that 
case, 

r(r) = M+Td. 

c 

The coefficient a/c(i^) hence read: 



^ + ^+i2^/c 
The related characteristic equation reads: 

1 _ e-^^^(l-e-(i + ^)) 

These equations are relatively complex to solve analytically in that general form (see e.g. [57]). 
However, when considering purely synaptic delay case (corresponding formally to c = 00, 
i.e. disregarding the transport phenomenon), a similar analysis as the one performed in 
the previous section allows computing in closed form the curves of Turing-Hopf instabilities. 



^The stabilization effect already mentioned in previous examples comes from the fact that noise tends to 
make the sigmoidal transform less sharp, which explicitly appears in our equations, since the slope of the 
effective sigmoidal transform reads gj \/{2'k{ \ ^g^V{))) which is always smaller than the slope of the sigmoid 
that intervenes in the deterministic case, gj \f^^ and which decreases as noise is increased. 
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(a) Bifurcations of the synchro- (b) a = and A = 0.1, ICl (c) cr = 1 and A = 0.1, ICl 

nized system 



mm 

(d) cr = and A = 0.1, IC2 (e) cr = and A = 0.2, IC2 




(f) cr = and A = 0.1 



Figure 4: Stabihzation by noise of the spatiahy homogeneous state (fiQ^vo). (a) shows the bifurcation 
diagram as a function of a and A of the fuhy synchronized state, (b): for initial conditions ICl given by (22) 
and smah noise levels, the system stabilizes on a non-spatially homogeneous state (mode k = 1) but for 
larger values of A (not shown) or a (c) the state (i2o,vo) is stabilized and attractive, (d): same thing for 
initial conditions IC2 given by (23) shows mode k = 2 excited, and as A is increased (e), this state looses 
stability in favor of the mode k = 1 before {/^o^vo) is stabilized, (f) shows the instability of the mode related 
to k = 4. 



Dynamical instabilities occur for parameters such that: 1/52^4^2^ 
us denote by ujk the quantity: 



l/s)2 



> p-. In that case, let 



0^ 



the instability arises for parameters such that there exists m G Z for which is valid the 
relationship: 

r| = — ( — arctan(6>a;fe) — arctan(27r/c5) + 27rm ). 

Let us now numerically explore the bifurcations of this system. We aim at exhibiting 
different solutions that are not spatially homogeneous, corresponding to different wavenum- 
bers, and to this end specify particular initial conditions. We have seen that for small 
values of the noise parameters and large values of the slope ^, the fully synchronized system 
presented two different stable equilibria that we denoted (jJ^i^vi) and (/i2,'^2)- If all initial 
conditions belong to the attraction bassin of the same stable fixed point, or are equal to zero, 
then the only mode to be excited corresponds to /c = 0, and the neural field stabilizes on a 
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constant mode. If we now partition F into different connected zones and set up the initial 
condition so that neurons belonging to different subsets of F belong to the attraction bassin 
of the two different stable fixed points, in a balanced manner, higher modes are excited. As 
examples, setting the initial condition to: 



/i(r,t) = < 



1 r G [0,0.25], t e [-r,0] 

-Ire [0.75,1], t e [-r,0] (22) 
overwise 



we excite a non-constant mode /c = 1 as illustrated in Fig. 4(b), and setting the initial 
condition to: 

'l r G [0,0.05] U [0.5,0.55], t G [-r, 0] 
/i(r, t) = {-1 re [0.25, 0.3] U [0.75, 0.8], t e [-r, 0] (23) 
overwise 

we excite the mode A: = 2, as illustrated in Fig. 4(d). In these plots, the abcissa represents 
the location on the space F and the ordinate corresponds to time. As noise is increased, 
this mode looses stability in favor of the mode k = 1 first, before this mode looses again 
stability in favor of the spatially homogeneous solution (/io? ^o) as expected from the analysis 
of the stability of that fixed point. Higher modes prove relatively unstable, and illustrated 
in Fig. 4(f) for an initial condition corresponding to a wavenumber k = 4. After a short 
transient, the system stabilizes on a stationary solution corresponding to the mode k = 1. 

For axonal delays, the characteristic equation obtained is very complicated to solve an- 
alytically. These are simplified when considering wizard hat connectivity kernels, as shown 
in [46] . We can however obtain complex closed form expressions of the Turing-Hopf bifur- 
cation curves using [57]. These are not made explicit here, and simulations performed using 
an adaptation of the the code provided by Venkov in [58] yield very similar phenomena (not 
shown) . 

The dynamic Turing bifurcations in the particular one-layer case treated, chosen for 
simplifying analytical exploration, presented a limited set of spatio-temporal behaviors: we 
observed that only stationary solutions are found, either spatially homogeneous or charac- 
terized by a few typical modes. In particular, no wave or oscillatory activity was observed, 
because of our particular choice of connectivity function which was chosen of constant sign 
(here, positive) for biological reasons. In order to go beyond these behaviors, we now turn 
to the numerical study of a more complex neural field system composed of two layers in 
order to include both excitation and inhibition phenomena. We will observe in these cases 
complex irregular dynamic Turing patterns are observed in a particular range of noise levels. 

4.2.2. Dynamic Turing patterns in a two layers neural field with periodic boundary conditions 
We now consider a two layers neural field composed of an excitatory (labeled by the index 
1) and an inhibitory (labeled by an index 2) populations distributed over the space F = 
or F = [0,1]. Both layer produce currents sent to both the excitatory and the inhibitory 
layers, and the interconnection strength depends on the (functional or anatomical) distance 
between themselves in F (see e.g. [59]). We consider here constant sign connectivities with 
an exponential shape, with a width depending on the population: Ja{r^r') = e~l^~^ '^^''/^a 
for a = 1 or 2. The choice of the typical spatial extension of the kernel depends on type of 
modeling chosen. For instance, if one considers that F is a functional space (i.e. each r G F 



Stochastic Neural Fields Dynamics 



31 



corresponds to a function of the population, for instance the preferred orientation of the 
column), then local interactions are dominated by excitation and the connection between 
two very distinct functions (e.g. orientations) are dominated by inhibition, which motivates 
a choice si > S2. If F models the spatial location of each neuron, it is known that inhibitory 
connections are generally characterized by shorter axons and excitatory synapses project 
their synapses further, which would motivate a choice of typical sizes such that S2 > si. 

In the present section, we consider periodic boundary condition. This choice hence is 
close from functional neural fields settings, and for instance classically models orientation 
preference. The neural field can be seen as defined on the torus S^. Two other types of 
connectivity will be dealt with in Appendix B: {i) reflective boundary conditions in which 
the solution is virtually evenly continued at the boundaries and 1 and the convolution is 
done on R instead of [0, 1] and (ii) zero boundary conditions where the convolution only 
occurs on [0, 1] (which would correspond to a convolution on R with a null continuation 
outside the interval [0, 1]). 

We define the type function of a neuron (excitatory or inhibitory) by the function 
iy{i) G {1,2}. The connectivity kernels J^y are multiplied by a typical connectivity coef- 
ficient between the different populations, Wj^j^^ for the deterministic interactions, chosen, 
consistently with the previous finite-population example, equal to: 

/ 15 -12 \ 
^=^16 -5 J- 

Noisy interaction ajy^j^'ir^r') are considered equal to a coefficient cFj^y multiplied by the 
exponential kernels J^/ (r, r'), and the typical time constants Ou{t) of all neurons is considered 
equal to 1. In a network composed of N^j, neurons of type v in the population located at 
r^ G r, the equation of neuron i of type v{i) = a G {1,2}, in population a at location Tq, 
reads: 



/ pm 2 \ 

dV\t) = i-V\t)^Ia{r^,t)^Yl T^JT S WauMra.r^) S{r^,V^ {t-r{r^,r^)))\ dt 

P{N) 2 \ 

+ E E ]^ E ^a.Ur^.r,) S{r,,V\t-T{r^,r,))) dB^^^^^^^ ^K{r^,t)dWl 

A straightforward multi-dimensional extension of theorem 2 and proposition 6 allows proving 
that the same propagation of chaos property applies, deriving the stochastic mean-field 
equations, proving that the solutions of these are Gaussian, with mean and covariances 
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satisfying the coupled equations: 

' ^(r, t) = -/XI (r, t) + \{r')dr'{wnJi{r, r')f{r, fi^ (r', t - T(r, r')), «i (r', t - T(r, r'))) 
+w'i2 J2(r, r')/(r, //2(r', t - T(r, r')), W2(r', t - T(r, r')))} + /i(r, t) 
^(r, t) = -Ai2(r, t) + A(r')dr'{«;2i Ji(r, r')/(r, ^ii{r', t - T(r, r')), vi{r' , t - T{r, r' j)) 
+W22 J2(r, r')f{r, ^2{r', t - T(r, r')), 1-2 (r', t - T(r, r')))} + /2(r, 
' If (r, = -2 «i(r, + /j. A(r')2dr'{<72i J2(r, r')p{r, t - T(r, r')), «i(r', t - T(r, r'))) 

+(t22 J|(r, r')/2(r, M2(r', t - T(r, r')), V2(r', t - T(r, r')))} + A2(r, t) 
If (r, t) = -2 V2{r, t) + A(r')2dr'{aii J2(r, r')f{r, ^i(r', t - T(r, r')), t;i(r', t - T(r, r'))) 
+cTi2 J|(r, r')/2(r, M2(r', t - r(r, r')), t;2(r', t - r(r, r')))} + Ai(r, t) 

(24) 

We numerically explore the solutions of these equations and the dependence upon noise, 
initial datum and boundary conditions on the neural field. We further simplify the problem 
by considering that the additive noise parameters Ai(r, t) and A2(r, t) are identical and 
constant, and denote their common value A. Similarly, we assume that cTj^^y/ are all equal to 
a quantity denoted a. The (space-dependent) delays are neglected in the present study, i.e. 
set T{r^r') =0. Their specific infiuence will be analyzed in a separate study. 

Since the connectivity matrix is convolutional and the domain is periodic, the system 
clearly satisfies the total synchronization condition given in equation (6), and we expect 
to find spatially homogeneous solutions when the initial conditions are themselves spatially 
homogeneous. Let us denote by Ja the quantity: 

Ja = j^Kr')dr'Ja{r,r') 

and by J^: 

for a G {1,2}. The infinite dimensional system (24) reduces, when considering its spatially 
homogeneous solutions, to the 4-dimensional ordinary differential equation: 

Vi = -111 +^iiJi/(/ii,vi) + ^12^2/(^2,^2) ^h{t) 

1^2 = -^2 +^2iJi/(M1,^i) +^22^2/(^2,^2) -^h{t) 

' V, =-2^i+a2(jf/2(/ii,^i) + j|/2(/i2,^2))+A? (2^) 

= -2^2+^2 (jf/'(/il,^l)+ J2'/'(/i2,^2)) +A2 

Functional Connectivity Case. We start by considering a case where Si > 52, which is 
related as stated to functional connectivity: local interactions are governed by excitation 
and distal interconnections by inhibition. We fix si = 0.02 and 52 = 0.0125. Figure 5 
represents the bifurcation diagram of the fully- synchronized system for /i = and I2 = —3, 
as a function of the additive noise amplitude Ai = A2 = A, and of the constant common 
synaptic noise amplitude a — a^yf. The codimension two bifurcation diagram presents a 
Hopf, a saddle- node and a saddle- homoclinic bifurcation curves, separating the diagram 
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(a) Codimension 2 diagram (b) cr = 0.1 (c) A = 0.1 

Figure 5: Bifurcation diagram of the fully synchronized state, (a) Codimension 2 bifurcation diagram with 
respect to a and A presents three bifurcations manifolds: saddle-homoclinic (SH, green), saddle-node (SN, 
blue) and Hopf (H, brown) separating the bifurcation diagram into 4 zones: two stationary (cyan, A and 
D), a periodic zone (C, orange) and a bistable zone (B, yellow). Codimension 1 bifurcation diagram as a 
function of A or cr for fixed values of the other noise parameter are displayed in (b) and (c). 



into three qualitatively distinct zones: two stationary zone (A and D, colored in cyan) 
where the moment system presents a stable fixed point, i.e. the solution of the mean- 
field equation is stationary, separated by the saddle-homoclinic bifurcation curve from a 
bistable zone (B, yellow) where the moment equations present both a stable fixed point 
and a periodic orbit, and a zone where there only exists a periodic orbit (C, orange). The 
bistable zone is separated from the periodic behavior zone through a saddle-node bifurcation, 
and the periodic behavior zone from the stationary zone (D) through a Hopf bifurcation. 
The analysis of this diagram underlines the similar but asymmetrical role of a and A. As 
examples we plotted codimension 1 bifurcation diagrams for A = 0.1 as a function of a and 
for <j = 0.1 as a function of A (black lines in the codimension two diagram). We observe 
that for small noise amplitude (either a or A) corresponding to parameter region (A), the 
system presents a stable high fixed point and two unstable fixed points, a saddle and a 
repulsive fixed point. As this noise parameter is further increased, a cycle appears from 
the saddle fixed point through a saddle homoclinic bifurcation, and we enter the bistable 
(yellow) region (B). As noise is further increased, the stable and the saddle fixed point merge 
and disappear through a saddle-node bifurcation, and the system is left with an unstable 
fixed point and a stable periodic orbit (region C). This orbit decreases amplitude disappears 
through a Hopf as noise is increased, and at this bifurcation the unstable fixed point gains 
stability (region D). 

The diagrams presented in Figure 5 characterize the existence and the nature of the syn- 
chronized states. As soon as the initial condition is homogeneous, the system will present 
solutions that are constant in space, and their time profile is given by the solutions of the 
fully-synchronized system. However, these synchronized states might not be stable, and in- 
homogeneities in the initial condition might lead the system to different states. Bifurcations 
from these spatially homogeneous states could be analytically as outlined in the previous 
section. However, due to the dimensionality of the system and to the particular coefficients 
chosen, the equilibria are complicated to characterize and hence the analytical characteri- 
zation of Turing-Hopf bifurcations is here intricate. This is why we address that problem 
numerically by computing the solutions of the neural-field moment equation for different 
values of the parameter A or a, with homogeneous or inhomogeneous initial condition. Re- 
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suits are displayed in Figures 6 and 7. In these figures is plotted the mean of the activity of 
the excitatory population, /ii(r, t), as a function of r (abscissa) and t (ordinate). 

The behaviors identified through the bifurcation analysis of the synchronized equations 
are perfectly recovered in the analysis of the neural field when taking an homogeneous 
initial condition, as expected from proposition 3. The stability of these synchronized states 
is analyzed through perturbations of the synchronized initial condition. Let us start by 
analyzing the effect of the parameter A for a fixed value of a = 0.1. In the stationary 
regime (A) corresponding to low noise amplitudes, very intricate phenomena take place that 
tend to destabilize the spatially homogeneous solutions. Let us fix for instance the initial 
condition of fii{r) to 5 for r G [0,0.05] and otherwise^. If the noise is small enough, 
the initial condition creates bi-directional waves that travel through the neural field. These 
waves split into different secondary waves, themselves potentially splitting. All these waves 
interact together, and this phenomenon results in highly irregular transient behaviors (see 
e.g. Fig. 6(a) for A = 0.1). This wave splitting and interaction phenomenon becomes 
sustained in time as A is further increase. For our parameters, the wave splitting and 
interaction phenomenon is sustained for A > 1. Figure 6(d) illustrates this behavior for 
A = 1. From the analysis of the spatio-temporal patterns, we observe the fact that waves 
tend to spread increasingly fast across the neural field as noise is increased. The irregular 
pattern suddenly turn as noise is further increased into a space-time quasi-periodic waves. 
These states, presented in Fig. 6(f), show more regular quasi-periodic spatial patterns that 
oscillate quasi-periodically in time. These waves appear not to interfere together, which 
explains the increased regularity observed in contrast with the patterns observed for smaller 
values of the noise parameter. These are present for increasing values of A until the bistable 
parameter region (B) is reached. 

For noise levels corresponding to the bistable region (B) or to the oscillatory region (C), it 
appears that the spatially homogeneous periodic regime is attractive: for non- homogeneous 
initial conditions, the system smoothly goes from a non-synchronized irregular activity to 
spatially homogeneous time-periodic activity corresponding to the spatially homogeneous 
solution of system (25) identified, see Figure 7(a). In the bistable region (B), one may wonder 
if similarly to the case of section 4.2.1, the presence of two stable equilibria could excite higher 
wavenumber. In order to numerically investigate this question we set up the initial condition 
at equilibrium for one part of the neural field and on the limit cycle for the other part. In the 
extensive numerical tests performed, we always obtained that the oscillatory regime would 
take over and govern the dynamics after very short transients, even if most of the neural 
field has its initial condition in the attraction basin of the spatially homogeneous stable fixed 
point. The unconditional stability of the spatially homogeneous behavior is also recovered 
for values of A in the stationary region (D): the spatially homogeneous stationary state is 
attractive for the system with inhomogeneous initial conditions and the system converges 
towards a constant spatially homogeneous solution (see Fig. 7(c)). 

For fixed values of A and a value of a varying, the exact same type of instability of the 
spatially homogeneous equilibrium in region (A) and stability in regions (B)-(D) appears 
(see Figure 8). Similarly, instabilities in region (A) is characterized by irregular wave inter- 
actions. The main difference between the effects of A or a lies in the fact that in the former 
case, the standard deviation of the solution of the mean-field equation has a constant value. 



^We observed in our numerical simulations that the qualitative phenomena described do not sensitively 
depend on the choice of the initial condition. 



Stochastic Neural Fields Dynamics 



35 





(c) A = 1, homogeneous IC (d) A = 1 (e) A = 1 



150 



200 
t 



(f) A = 1.6 



Figure 6: Irregular transient phases for inhomogeneous initial conditions, for parameters in region (A): 
fii(r,t) solutions of the two-layers mean-field equations is plotted. Spatio-temporal patterns: Abscissa: 
location on F = [0,1]. Ordinate: time. Plots (b) and (e) represent /ii(r = 0.1, t) (black) and /ii(r = 0.5, t) 
(red). Initial conditions are set to /ii(r, 0) = 5 for r G [0, 0.05] and elsewhere (orange box) /i2(^, 0) = and 
v{r,0) = 0, except in (c) where we show the solution for spatially homogeneous (null) initial conditions, (f): 
Space-time oscillation: A = 1.55. Waves no more interfere and the activity forms a non-stationary quasi- 
periodic structure. Right: plots of the activity for two different times as a function of r G F. Animations 
of the activity are available in the supplementary material. Figures and animations were obtained using 
XPPAut [60]. 
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Figure 7: Spatially homogeneous stable behaviors: A belongs to the bistable (B), periodic (C) and stationary 
(D) regimes. Same inhomogeneous initial conditions and setting as in Fig. 6. Animations of the activity are 
available in the supplementary material. 
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(d) cr = 1.6, Variance (e) Different trajectories (f) a = Q 

Figure 8: Behaviors of the two-layers mean-field equations for different initial conditions and synaptic noise 
amplitude a. (e): 100 trajectories from the fully-synchronized equations stochastic equations in a network 
of 10 000 neurons show a sharp synchronization when the standard deviation reaches zero (blue: excitatory 
neurons (population 1) and red: inhibitory neurons (2)). 

i.e. each individual trajectory in the network deviates from the mean activity by a Gaus- 
sian process with constant standard deviation, whereas in the latter case corresponds the 
standard deviation of the solution of the mean-field equations depends on the mean activity 
of the neural populations. This phenomenon is particularly visible in the non-stationary 
solutions. For instance for a = 1.6, we observe that the standard deviation reaches very 
small values at times close to the onset of an oscillation, as shown in Figure 8(d), resulting 
in the fact that all neurons sharply synchronize at these times, a phenomenon particularly 
visible in the network simulations plotted in Fig. 8(e) where we plotted 100 trajectories of 
a network with 10 000 neurons. 

All this analysis illustrate the fact that noise strongly shapes the form of the activity 
of large network, governing the nature of spatially homogeneous states, and their stability, 
yielding different spatio-temporal patterns that can be either transient or sustained. Irreg- 
ular behaviors arose exclusively for relatively small values of the noise parameter, which 
again illustrates a stabilizing effect of noise in favor of the spatially homogeneous in law 
solutions. Very interesting irregular patterns of activity characterized by the irregular for- 
mation of waves and their splitting was observed in parameter regions close the presence of 
a saddle-homoclinic bifurcation in the fully synchronized system. This same phenomenon 
is found for similar systems defined on [0, 1] with different boundary conditions: reflective 
or zero. Results are displayed in Appendix B. This formation of complex spatio-temporal 
patterns arising from a wave-splitting phenomenon strongly evokes Turing patterns as found 
in different reaction diffusion equations in biological mathematics. In particular, these pat- 
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terns are similar to those exhibited by [61] and obtained from the analysis of the dynamics 
of reaction diffusion equations related to pattern formation on the shells of mollusks, in a 
case where the system induces the formation of forward and backward running waves that 
interact. These can be analyzed as done in [62] in a case where are identified self-replicating 
bumps, compared with dynamics observed dissipative equations such as Ginzburg-Landau's, 
but to our knowledge not observed in neural field equations. This structuring effect of noise 
and the presence of the particular type of spatio-temporal patterns observed in that system, 
beyond illustrating the fact that noise levels strongly shape the dynamics of neural fields. 

Anatomical Connectivities: Bumps, bump- splitting and interfering waves. In this section we 
consider the anatomical case where distal interconnections are dominated by excitation and 
local interactions by inhibition. In that case, the codimension one bifurcation diagram of 
the synchronized system is provided in Figure 9(a) and shows the same global structure. 
This diagram characterizes the nature of spatially homogeneous solutions. We investigate 
solutions of the mean-field equations corresponding to non-homogeneous initial conditions, 
which, as observed in the previous sections, do not necessarily converge towards spatially 
homogeneous patterns. 

The observations are essentially of the same kind as identified in the functional connec- 
tivity case. First of all, the bifurcation diagram of the spatially homogeneous equations 
can be segmented, like the previous case, into the same qualitative regions (A)-(D), and it 
appeared in that case again that the spatially homogeneous states were unstable only in 
the small noise region (A). However, we observed that depending on the ratio r = Si/s2, 
different kind of spatio-temporal patterns appeared (see Figure 9). For ratios greater than 
r* ~ 0.6, the same phenomenon of wave splitting and interactions as in the functional case 
is observed (see Fig. 9(1)). 

For r < r*, a new type of dynamics appear, characterized by the presence of bumps of 
activity, i.e. localized stationary patterns. These are evidenced for initial conditions equal to 
zero except on an interval [a, b] G S^. For A = 1, non spatially homogenous solutions appear, 
composed of a specific number of bumps that depends on the width of the interval [a, b]. As 
noise is increased, these bumps tend to split in two different bumps of the same spatial size 
and either stabilize, or split again, a phenomenon strongly evocative of the patterns observed 
in a different context by Coombes and Owen in [62] . This sequence of bump splitting can 
either stabilize into a stationary pattern composed of several bumps or appear repetitively 
depending on the value of the noise parameter. For relatively small noise values (A < 1.7), 
the activity stabilizes on a stationary pattern which is not spatially homogeneous, and the 
number of bumps and their shape depend both on the initial condition chosen and on the 
value of the noise. It is interesting to note that we observed that for A > 1.5, the stationary 
pattern found is spatially periodic, characterized by a specific wavenumber increasing as 
noise is increased, and depending on the type of initial condition chosen. For instance, we 
show in Figure 9(e) the case of A = 1.5 and initial condition zero except on [0,0.05] where 
it is equal to 5 (ICl), where the stationary behavior is characterized a spatially periodic 
pattern with wavenumber 8, and for initial condition zero except on [0,0.05] U [0.1,0.15] 
where it is equal to 5, the same kind of phenomenon appears and stabilizes on a pattern 
with wavenumber equal to 9. For A = 1.6 and initial conditions ICl, the wave number is 
10. For values of A larger than 1.7, the system no more stabilizes into stationary patterns 
and presents very complex irregular sequences of bump splitting, before presenting irregular 
spatio-temporal waves (A = 2.2, Fig. 9(i)). 
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(a) Spatially Homogeneous 



(b) A = 1 



(c) A = 1.2 




(d) A = 1.3 



(e) A = 1.5 



(f) A = 1.5,IC2 





(g) A = 1.6 



(h) A = 1.7 



(i) A = 2.2 






(j) A = 2.3 



(k) A = 2.5 



(1) A = 1.6, ^ > r* 



Figure 9: Anatomical connectivity case with /ii(r, t) represented as a function of r G S-*^ (abscissa) and t 
(ordinate), (a)-(k): si/s2 = 0.57 show a sequence of bump splitting as A is increased in the parameter 
region where stationary spatially homogeneous solutions related to small values of A. Orange Box: initial 
condition /ii(r, 0) = 5 for r G [0,0.05] and otherwise, and blue box: — /ii(r, 0). (1): si/s2 = 0.65: wave 
splitting phenomenon, for A = 
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These bumps disappear in favor of a spatially homogeneous periodic activity when noise 
levels reach the bistable region, and this region turns into spatially homogeneous stationary 
solutions as noise is further increased, in the case corresponding to regions (B), (C) and (D). 

For ratios Si/s2 greater than r*, as stated, s we find phenomena resembling to what was 
found the functional connectivity case: for small values of the noise parameter, a wave- 
splitting and interference phenomenon is found, progressively turning into quasi-periodic 
spatio-temporal waves, spatially- homogeneous oscillations and spatially homogeneous sta- 
tionary solutions as the noise is increased. The wave splitting phenomenon is displayed on 
Figure 9(1), other cases resemble to the other patterns found and are not displayed. 

5. Discussion 

In this article, we addressed the problem of characterizing the dynamics of large as- 
semblies of neurons that are spatially extended and interact after delays. Modeling each 
individual neuron and characterizing their interactions, we derived neural field equations 
using a rigorous probabilistic mean-field theorem developed in [36]. The equations that 
were obtained in that article were implicit equations on stochastic processes probability dis- 
tributions, and are extremely hard to tame in that general framework. We summarized the 
mathematical results obtained in a general setting, and provided a simple condition for the 
existence of solutions that are spatially homogeneous in law. The main result of the present 
article is the precise analysis of the dynamics of these equations in the particular case of 
firing-rate neurons, corresponding to the fact that membrane potential of the neurons have 
a linear intrinsic dynamics with nonlinear interactions. We were able in that case to re- 
duce the mean-field equations to a set of delayed differential or integro-differential equations 
that are mathematically tractable, and using this description evidenced the important role 
of the microscopic noise in the dynamics in a series of examples of increasing complexity: 
finite-populations networks with delays, one-layer neural fields and two-layers neural fields. 
Several phenomena such as cascades of Hopf bifurcations, Turing-Hopf bifurcations and ir- 
regular spatio-temporal patterns were evidenced for specific levels of noise, and different 
levels of noise were identified corresponding to specific typical behaviors. 

The first obvious limitation of the study is the fact that this precise analysis reduces to 
the case of firing-rate neurons. Though popular in the study of neural fields in the neu- 
roscience community and widely used, this choice can be discussible in the present study. 
Indeed, the model does not takes into account the highly nonlinear nature of several neu- 
ronal phenomena. This linear approximation is often seen as a heuristic mean-field limit 
itself, which can raise questions of the relevance of our present approach for neuroscience ap- 
plications. The analysis developed here might hence better account for large-scale neuronal 
areas in which a very large number of cortical columns, described by the Wilson and Cowan 
model, s interact. As far as neurons are concerned, we may add that a large body of experi- 
mental results have proved successfully accounted by linear dynamics, generally relative to 
the specific nonlinear ly determined fixed point of normal brain activity. These linear models, 
though less general and accurate that the nonlinear ones, yielded much greater insight, and 
in particular analytic treatment, of the cortical dynamics within their regimes of validity. 
Probably the main reason why we chose this type of models is the reduction allowed by the 
Gaussian nature of the solutions. Indeed, as we demonstrated here, these models have the 
great advantage to be exactly and rigorously reducible to a set of tractable deterministic 
equations. This is hence a breach to go further in analyzing the complex dynamics of the 
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neural fields mean-field equations, with the aim of further understanding nonlinear neuron 
models. The analysis provided here is the first (and to our understanding, the only case 
where such an analytical study is possible) to address precisely the dynamics of such com- 
plex mean-field equations, and can also be seen as a proof of concept of the dynamics of this 
class of equations, in particular the effects of noise and delays in these equations. All these 
reasons motivated the investigation of linear models to uncover the dynamics of the more 
complex nonlinear neural fields models. We expect that nonlinear equations of this kind will 
also reveal other behaviors that are not reproduced in our linear approach. A perspective 
of great interest also would be to derive from the non-linear mean-field equations systems 
governing macroscopic variables such as the mean firing-rate. This is a complex and deep 
question we are currently investigating. 

The general version of the results of theorems 1 and 2 have several implications in 
neuroscience that are discussed in [36]. Among these, we may recall that the propagation 
of chaos property implies an independence between the behavior of finite sets of neurons, a 
surprising property in the classical view of neural populations postulating that the activity 
of neurons belonging to the same neuronal population, since sharing common input and 
strongly interconnected, was highly correlated. This view was recently contradicted by 
several experimental articles that exhibited high quality data corresponding to extremely 
low levels of correlation, in the primary visual area VI of awake macaques [10] and in the 
rodent neocortex [11]. These results confirm that the chaotic state is biologically relevant, 
and also have several implications in neural coding, as shown in the analysis of the coding 
efficiency of decorrelated activity in [10]. The independence of the neural activity also 
reconciles spike [9, 63] and rate [64] coding, two opposing conceptions of the neural activity, 
by emphasizing the role of collectivity: observing the state of several independent neurons 
during a time corresponding of the emission of a few spikes allows precise evaluation of the 
probability distribution of the interspike interval, hence the firing rate of neurons. 

One of the main new feature of the mean-field model is the consubstantial integration 
of the noise in the dynamics of the neural field: in the reduced equations obtained for the 
linear systems, we observed that noise appeared as a parameter in the sigmoidal function 
transforming the voltage into a mean- firing rate. This property quantified precisely the ef- 
fect of noise in the behavior of neural fields through our mean-field equations. We studied 
a finite-population mean-field model with delays, and showed noise and delays shape the 
behavior of the system. As delays are increased, a sequence of Hopf bifurcation appears in 
the system, that was observed numerically and accounted for analytically by the analysis 
of the characteristic roots for a particular system. This sequence of Hopf bifurcations pro- 
duces complex transient behaviors, and simultaneously with noise levels shape the nature 
of the response of the neural field, by producing either stationary or oscillatory responses. 
Continuous neural fields were then analyzed with the particular aim of exhibiting qualita- 
tive changes in the behavior of a neural field due to noise levels. This is why we addressed 
the presence of Turing-Hopf bifurcations, first in a one-layer system semi- analytically, and 
then numerically in a more complex two-layers neural field. We evidenced several noise- 
induced dynamical Turing-Hopf bifurcations generating different spatio-temporal patterns, 
sometimes irregular and extremely complex. We also observed that noise also governs the 
stability of synchronized in law behaviors. 

The periodic solutions triggered by an increase in the noise standard deviation corre- 
spond at the level of the network to solutions in which in neuron fires synchronously, in 
systems that present a unique stable fixed point when neurons are not driven by noise. 
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Noise hence appear to have a very strong structuring effect on the dynamics of the sys- 
tem, which can appear at first sight counterintuitive as noise is usuahy seen as altering the 
structure of the activity. These observations add to different phenomena documented in the 
neural networks literatures, as for instance coherent oscillations in neural networks [65, 66] 
where evidences for the presence of noise-induced in noisy integrate- and- fire neurons. The 
phenomena observed in our analysis differs from the more classical stochastic resonance or 
coherence resonance phenomena well documented in the neurocomputational literature (see 
e.g. [67] for a review of the effect of noise in excitable systems), since these correspond to the 
fact that there exists a particular level of noise maximizing the regularity of an oscillatory 
output related to periodic forcing (stochastic resonance) or to intrinsic oscillations (coher- 
ence resonance). Another important result that can be derived from the present analysis 
is the ability to define classes of parameter ranges attached to a few generic bifurcation 
diagrams as functions of the parameters. This property suggests model-based strategies for 
parameter evaluations in experimental measurements. 

The present analysis, reducing noise levels to a simple parameter of a dynamical system, 
allowed going further into the analysis of the functional role of noise in the brain, a question 
currently widely debated in the neuroscience community, and relates this question to the 
presence of synchronized oscillations, as we showed in our different numerical simulations 
that noise levels determined the presence of oscillations, and the type of noise (additive 
or synaptic) the sharpness of the synchronization of oscillations across different neurons or 
areas. 

Several extensions of the present study with applications in neuroscience and in applied 
mathematics are envisioned. The first would be, as already mentioned, to investigate similar 
phenomena in biologically relevant models and to study the behavior of the solutions of the 
mean-field equations in that mathematically much more complex setting that would include 
in particular nonlinear, excitable intrinsic dynamics and different ionic populations. The 
mean-field equations in that setting are complicated implicit stochastic equations, that can 
be written as a non-local nonlinear partial differential equation on the space of probabil- 
ity distributions. These equations are extremely hard to analyze directly, and the proof of 
concepts provided by this simpler analysis foreshadows interesting behaviors. An additional 
intricacy of the model is that, as we observed in our analysis, the system does not always 
present stationary solutions, and when these exist, they are not necessarily stable. Cycles 
were found, and bistable states between stationary and periodic behaviors were found, as 
well as irregular states in spatially extended systems. Analyzing the solutions of the non- 
linear mean-field equations will hence present interesting technical difficulties and probably 
complex behaviors, and necessitates the development of techniques to analyze the conver- 
gence of such stochastic equations towards periodic solutions. Another important direction 
in the development of this work would consist in fitting the microscopic model to biological 
measurements. This would yield a new neural assembly model for large scale areas and de- 
velop studies on the appearance of stochastic seizures and rhythmic activity in relationship 
with different parameters of the model, integrating the presence of noise in a mathematically 
and biologically relevant manner. 

A question that remains largely open is the question of the number of neurons per pop- 
ulation. Motivated by the fact that typical sizes of the neuronal populations are generally 
orders of magnitude larger than the number of neural populations (see e.g. [44]), we followed 
the choice made in [36] and assumed that the number of neurons in each population was 
increasing fast enough with the number of neurons (assumption of equation (3)), which left 
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room for the propagation of chaos to occur in each population. If this condition is not sat- 
isfied, different hmits can appear, and characterizing these different hmits is a chahenging 
mathematical question. Another very interesting and deep question in this setting is the 
question of mean-field limits in slow-fast systems, problem that is highly relevant for differ- 
ent aspects of the neuronal behavior, and which raises a number of technical and applied 
questions. 



Appendix A. Existence and Uniqueness of solutions of the synchronized mean- 
field solution 

In this appendix we prove the existence and uniqueness of solutions of the fully synchro- 
nized system given by equations (7). This proof is similar to the proof given in [36]. Since 
the quantities 5(r, x, cp) and i^(r, x, (f) do not depend on r, we drop the dependence of these 
functions in r. As usually done, we transform the equation (7) into a fixed point equation 
on the space of stochastic processes. To this end, let us define the map <l> as follow: 



\S + J^ds^G{s,Xs)^Mz[B{Xs,Z^.^)] 
[{Zt) ={Xt) e M independent of (X,), {Wt{-)) and •)) 



The solutions of equation (7) are exactly the fixed points of ^. We assume here that G 
and g are i^-Lipschitz-continuous and satisfy the linear growth condition, and b and (3 are 
L-Lipschitz continuous in both their variables. It is easy to show that any possible solution 
has a bounded second moment following [36] . 
Existence- 
Let X^ e M^{C) such that XO|[_^o]=Co a given stochastic process. We introduce the 
sequence of probability distributions (X^)/e>o on A4{C) defined by induction as = 
(<I>(X^)). We denote by (Z^) a sequence of processes independent of the collection of 
processes (X^) and having the same law. We decompose into six elementary terms the 
difference: 



X, 



= ^ (g(s, X^) - G{s, Xt')) ds 

B{X^,,Z^)- B{X^,-\Z^) 



E; 



ds 



B{X^-\Z'') - B{X^^-\Z^-^) 



f 

Jo 

^ [g{s,X^)-g{s,X^-'))dW, 
H{X^,Z'')-H{X^-\Z'' 



ds 



■ E^ 



H{X^-'^,Z'') - H{X^-^,Z''-'^) 



def 



At + Bt + Ct + Dt + Et + Ft 
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where we simply identify each of the six terms A^, ^t, C^, Dt, Et and Ft with the corre- 
sponding expression in the previous formulation. By a simple convexity inequality (Holder) 
we have: 

and treat each term separately. 

The term At is easily controlled using Cauchy-Schwarz inequality, Fubini identity and 
standard inequalities and we obtain: 









E 


sup \As\^ 








Jo ^ 



sup \Xt-Xt-'\^ 

-r<n<s 



ds 



Similarly, the martingale term Dt is bounded using the Burkholder-Davis-Gundy theorem 
to the d-dimensional martingale ( Jq {g{s, X^) — g{s^ dWs) and we obtain: 



E 



sup \Ds 

0<s<t 



<4.K^ / E 



Jo 



sup \x!:-xt 



k-l\2 



ds 



Let us now deal with the deterministic interaction terms Bt and (7+. We have: 



\Bt\' 



I ds I \{r')dr' f drjir, r', u)OEz[bir, r', X^ J - b{r, r', Xt\ J]) 
Jo Jr J-T 

Jo Jr J-T 

< tXirfK^L"^ f \X^, - |2 ds < tXiTfn'^L'^ f sup \Xl - Xl-^\^ ds 

Jo Jo —T<u<s 

hence easily conclude that 

^sn^ \Bs\^]<t\{TfK^L^ [\[ sup \X^-Xl-^\^]ds 
se[o,t] Jo -T<u<s 

and similarly 

E[sup <a(r)2/,2L2 ^Ei sup 

se[o,t] JO -T<u<s 

Eventually, the terms Et and Ft use Burkholder-David-Gundy (BDG) inequality in place 
of Cauchy-Schwarz' together with similar arguments as used for Bt and Ct. Using the 
cylindrical nature of the Brownian motions {Bt{r^r'))^ BDG inequality yields, for the term 
Et (and similarly for the term Ft): 

E[sup |e,p] <4A(r)2/^2^2 /'e[ sup \xl-x^-^\'']ds 
se[o,t] Jo -T<u<s 

for Qt equal to Et or Ft. Putting all these evaluations together, we get: 



E 



sup \X^^^-X^\^] <6{T^4){K^^2X(r)^L^K^) [ E[ sup \X^-X^-^\^]ds (A.l) 
se[o,t] -I JO -T<u<s 



Stochastic Neural Fields Dynamics 



45 



Moreover, since 



/c+l _ 



for t G [— r, 0] by definition, we have, noting 



=E 



sup |X; 

-T<S<t 



the recursive inequahty M/^ < M^'^ ds with = 6{T ^4){K^ ^2X{r)^ ^ i^^), which 

classicahy ahows concluding on the existence and uniqueness of solutions. 



Appendix B. Spatio-temporal patterns for one-dimensional neural fields with 
reflective or zero boundary conditions 

In this appendix, we reproduce the results obtained in section 4.2.2 with different bound- 
ary conditions. The system considered is essentially the same, i.e. a two layers excitatory- 
inhibitory neural field with exponential connectivity kernels, and set the connectivity matrix 
as done in the main text. In section 4.2.2, the neural field was defined on S^, i.e. character- 
ized with convolutional interactions and periodic boundary conditions on the connectivity 
kernel. We now choose in the present section different boundary conditions: reflective^ in 
which case the total synchronization conditions are satisfied, and zero boundary conditions, 
in which case no spatially homogeneous solution exists. 

Appendix B.l. Reflective boundary conditions 

In the case of a connectivity defined by an exponential connectivity kernel on F = [0, 1] 
with reflective boundary conditions, the stationary solutions are the same as in 4.2.2 and 
their bifurcation diagram given in Fig. 5. Numerical simulation of the spatio-temporal 
activity of this neural field for non-homogeneous initial conditions and different values of the 
noise parameter A are displayed in Figure B.IO. For small values of the noise parameter in the 
parameter region (A), the inhomogeneous initial conditions produce large amplitude waves 
that interact together, creating a transient complex structure of spatio-temporal activity 
that stabilizes on a fully synchronized stationary solution (i.e. constant in space and time). 
For slightly larger noise amplitude, this transient irregular phase takes over and produce 
sustained dynamic irregular activity on the neural field. As noise is further increased, a 
spatio-temporal periodic activity arises as the waves become faster and stop interacting. 
When the noise parameter reaches parameter regions (B) or (C), the spatially homogeneous 
oscillatory activity is recovered after a short transient. For noise levels in the parameter 
region (D), the whole neural field integrates fast the inhomogeneous initial condition and 
converges towards a spatially stationary activity. 

Appendix B.2. Zero boundary conditions 

In the case of zero boundary conditions of the connectivity kernel, the system does not 
fulfills the synchronization condition on the kernel given by equation (6). We analyze here 
the dynamics of this system for homogeneous initial conditions, and simulation results are 
displayed in Figure B.ll. In that case, neurons at locations r close to the boundaries and 
1 receive less input that neurons far from the boundaries, and non-synchronized patterns 
are observed for any parameter values. We however make similar qualitatively observations 
on the activity of the neural field. In details, for small noise intensities, the mean of the 
state of each neuron converge towards a constant value depending on its space location 
(see Fig. 11(a) for A = 0.1). As noise is increased, irregular spatio-temporal patterns of 
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(d) A = 1.5 (e) A = 3 



Figure B.IO: Solutions of the mean-field equations for the reflective exponential connectivity kernel, as 
noise amplitude is increased. Abscissa represents the location on F = [O?!]? Ordinate the time. Initial 
conditions are 5 on [0, 0.025] (orange box) and else equal 0. Animations of the activity are available in the 
supplementary material. 
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(c) A = 2 (d) A = 3 



Figure B.ll: Solutions of the mean- field equations for the exponential convolutional kernel with no (zero) 
boundary condition, as noise amplitude is increased. Abscissa represents the location on F = [0, 1], Ordinate 
the time. Initial conditions are deterministic equal to 0. Animations of the activity are available in the 
supplementary material. 

activity similar to these obtained in the previous cases are observed, and as noise increases 
this activity converges towards a periodic in time, non- synchronized activity that disappears 
into a constant in time, non-synchronized activity for larger values of the noise intensity. In 
this particular case, no transient irregular activity, nor spatio-temporal waves are found. 
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